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ABSTRACT 

An  analytical  approach  to  the  solution  of  the  whipping  response  of  a  submerged  submarine 
to  the  pulsation  of  a  nearby  explosion  bubble  is  presented    Particular  aspects  of  energy 
dissipation  (damping)  are  addressed  in  regards  to  their  effects  on  the  overall  loading  forces 
and  dynamic  response  of  the  submarine     Particular  damping  mechanisms  discussed 
include  hull  damping  (hull  material  and  structural),  internal  damping  (internal  structures 
and  sloshing  liquids),  and  external  damping  (hydrodynamic  damping  including  wave 
radiation  and  viscous  fluid  effects) 

A  flexible  analytic  model  is  developed,  using  the  fundamentals  of  vibration  and  hydro- 
mechanics.  A  finite-element  beam  model  is  used  to  represent  the  flexural  structure  of  the 
hull    Internal  structures  are  modeled  as  separate  dynamic  systems,  using  both  discrete  and 
modal  superposition  approaches    Onboard  liquids  are  modeled  using  a  mechanical- 
equivalent  system  based  upon  a  potential  flow  solution  of  the  liquid  free-surface    External 
hydrodynamic  forces  are  modeled  using  a  modified  Morison  formulation,  with  fluid 
velocities  and  accelerations  calculated  based  upon  a  popular  explosion  bubble  model    The 
equations  of  motion  for  the  composite  dynamic  systems  are  solved  in  the  time-domain 
using  a  modified  Newmark  direct  time  integration  scheme,  with  iterations  at  each  time 
step  accomplished  using  a  modified  Newton-Raphson  method 

The  analytic  model  is  implemented  on  a  personal  computer,  and  used  to  numerically 
investigate  various  damping  mechanisms  associated  with  dissipation  of  whipping  energy 
Additionally,  the  analytic  model  is  compared  to  explosion-induced  whipping  tests 
conducted  on  the  U.S.  Navy  test  platform  "Red  Snapper".  The  analytic  model  compared 
well  to  the  experimental  data  for  early-time  response  (through  the  2nd  bubble  period),  but 
tended  to  over  predict  response  for  later  times. 

Potential  weaknesses  with  the  analytic  model  are  discussed  in  light  of  the  comparisons 
with  experimental  data,  and  recommendations  for  future  analytic  work  in  the  area  are 
made 
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NOMENCLATURE 

V  vector  gradient  operator 

V  volume 

0  matrix  of  mode  shape  vectors,  velocity  potential  function 
a  Rayleigh  coefficient,  non-dimensional  bubble  radius 

(3  Rayleigh  coefficient,  viscous-frequency  parameter,  coefficient  for 

Newmark  algorithm,  control  variable  for  free  surface  effect 

5  logarithmic  decrement 

s  strain,  control  variable  for  bubble  migration 

y  coefficient  for  Newmark  algorithm,  adiabatic  constant  for  explosive 

r|  free  surface  elevation 

<J)  mode  shape  vector,  velocity  potential 

X  non-dimensional  rate  of  change  of  bubble  depth 

u  coefficient  of  friction 

v  kinematic  viscosity 

9  rotation  angle,  angular  orientation 

p  density 

a  stress,  non-dimensional  rate  of  change  of  bubble  radius 

x  non-dimensional  time 

co  radial  frequency 

Q  non-dimensional  bubble  depth 

C,n  modal  viscous  damping  ratio 

A  amplitude,  area 

C  damping  matrix 

CA  added  mass  coefficient 

CD  drag  coefficient  for  explosion  bubble 

Cd ,  Cj  drag  coefficient 

C  inertia  coefficient 

Cn  modal  damping 

D  specific  damping  energy,  diameter 

D0  initial  head  on  bubble 

E  Young's  modulus 

F  force  vector 

G  Green  function,  shear  modulus 

1  2nd  moment  of  area,  moment  of  inertia 
J  polar  moment  of  inertia 

K  stiffness  matrix 

KC  Keulegan-Carpenter  number 

Kn  modal  stiffness 

L  length 

M  mass  matrix 


M, 

mass  of  liquid  in  a  tank 

Mn 

modal  mass 

N 

number  of  cycles,  number  of  modes 

P 

pressure 

Qn 

modal  force 

Rn 

peak  Reynolds'  number 

S 

surface 

T 

period 

U0 

amplitude  of  relative  velocity 

V 

velocity 

w 

energy/work,  explosive  charge  weight 

x,x,x 

displacement  vector 

a 

tank  width,  bubble  radius 

c 

damping 

d 

bubble  depth 

e,,e. 

strengths  of  sources,  dipoles 

f 

force 

g 

acceleration  due  to  gravity 

h 

hysteresis  damping  constant,  tank  depth 

k 

stiffness,  constant  for  explosive 

m 

mass,  applied  moment 

n 

mode  number,  normal  vector 

q 

vector  of  modal  coordinates 

r 

radial  distance 

t 

time 

u,ii 

fluid  velocity,  acceleration 

V 

velocity 

Vm 

bubble  migration  velocity 

x,x,x 

displacement,  velocity,  acceleration 

1.   INTRODUCTION 

1.1        Underwater  explosions  and  ship  whipping. 

Designers  of  naval  ships  have  long  been  concerned  with  the  effects  of  underwater 
explosions.  Historically,  most  attention  has  been  paid  to  the  effects  of  contact  charges  and 
charges  of  near  to  intermediate  standoff   Protection  methods  to  these  effects  have 
typically  included  heavy  side-protective  systems  to  minimize  possibility  of  hull  rupture  and 
shock-mounting  of  vital  machinery  and  equipment  to  minimize  damage  caused  by  the  high 
frequency  shock  wave  resulting  from  the  explosion 

While  hull  rupture  and  internal  shock  wave  effects  are  still  considered  of  major 
importance  to  ship  survivability,  increasing  attention  has  been  paid  in  recent  years  to  the 
whipping  of  ship  hulls  caused  by  charges  of  near  to  intermediate  standoff  Whipping  can 
be  defined  as  the  transient,  flexing  motion  of  the  whole  hull  of  a  ship  in  its  vibrational 
modes  corresponding  to  the  lowest  natural  frequencies1.  Such  motion,  if  of  sufficient 
magnitude,  can  cause  hull  girder  buckling,  tearing,  or  other  loss  of  hull  girder  strength. 

Numerous  examples  of  catastrophic  or  near-catastrophic  whipping  can  be  cited 
throughout  naval  history.   As  recently  as  1991,  the  U.S.  Navy  Guided  Missile  Cruiser 
USS  Princeton  (CG-59)  struck  a  floating  contact  mine  in  the  Persian  Gulf.  Although  the 
local  hull  rupture  caused  by  the  mine  explosion  (near  the  bow)  was  not  a  threat  to  the 
survival  of  the  ship,  the  subsequent  whipping  of  the  hull  caused  near-catastrophic  hull 
girder  damage  near  the  stern  of  the  ship. 

When  an  underwater  explosive  detonates,  the  solid  explosive  material  suddenly 
reacts,  leaving  behind  gaseous  products  at  very  high  temperature  and  pressure    Almost 
immediately,  the  shock  wave,  produced  due  to  the  sudden  discontinuity  in  pressure, 
propagates  radially,  approximately  with  the  speed  of  sound  in  water  This  shock  wave 
carries  with  it  approximately  one-third  of  the  energy  released  during  the  reaction2    It  is 
this  shock  wave  that  is  typically  responsible  for  localized  hull  rupture  and  damage  to 
internal  equipment. 

Left  behind  as  the  shock  wave  propagates  through  the  fluid  medium  is  a  cavity  of 
gaseous  reaction  products  at  high  pressure.  This  cavity  subsequently  expands  (as  a 
bubble)  to  relieve  the  difference  in  pressure,  accelerating  the  surrounding  fluid    The 
bubble  continues  to  expand  beyond  the  point  of  hydrostatic  equilibrium  (due  to  the  inertia 
of  the  surrounding  fluid)  until  a  point  of  dynamic  equilibrium  is  reached.  The  bubble  then 


•Hicks.  A.N..  "Explosion  Induced  Hull  Whipping",  Advances  in  Marine  Structures.  Admiralty  Research 
Establishment.  Dunfermline.  Scotland.  UK.  1986.  p.  390. 
2Ibid.,p.  392. 


reverses,  continuing  to  contract  until  dynamic  equilibrium  is  again  reached,  where  it 
quickly  rebounds  and  again  begins  to  expand    This  bubble  expansion  and  contraction 
continues  until  the  energy  of  the  reaction  is  fully  dissipated     As  the  bubble  rebounds,  it 
greatly  accelerates  the  surrounding  water,  generating  a  substantial  pressure  pulse 
(commonly  known  as  the  bubble  pulse).    This  bubble  pulse  can  impart  significant  loads  on 
structures  in  the  vicinity 

Thus,  for  certain  underwater  explosions,  a  ship  can  experience  loading,  not  only 
from  the  radiating  shock  wave,  but  also  from  the  quasi-periodic  bubble  pulse.  Although 
the  whipping  motion  of  a  ship  hull  is  characterized  by  free  vibration  of  the  hull  in  the 
water,  it  is  this  quasi-periodic  loading  which  can  reinforce  the  free  vibration  and  magnify 
the  response  of  the  hull  to  whipping,  particularly  when  the  period  of  the  bubble  pulse  is  in 
the  vicinity  of  the  natural  periods  of  hull  flexure. 

1.2        The  importance  of  damping  in  ship  whipping. 

The  first,  most  easily  observed  effect  of  damping  on  structural  vibration  is  in  the 
decrease  in  amplitude  during  free  vibration    Figure  1-1   illustrates  the  effect  of  damping 
on  the  free  vibration  of  a  cantilever  beam.   If  the  cantilever  with  some  damping  is 
deformed  and  then  released,  it  will  oscillate  with  a  regular  period,  but  the  amplitude  of  the 
oscillation  will  decrease  with  time.  Without  damping,  the  cantilever  would  continue  to 
vibrate  without  the  decrease  in  amplitude.  In  fact,  the  free  vibrations  of  damped  structures 
will  normally  be  reduced  at  a  rate  that  may  be  used  as  a  measure  of  the  amount  of 
damping  in  the  system    A  well  known  measure  of  damping,  known  as  the  logarithmic 
decrement,  5,  is  related  to  the  ratio  of  the  nth  to  the  (//  -  A^th  cycle  amplitudes  by: 

5  =  -J-ln(-^-)  (1-1) 

N       An+N 

where  An  is  the  amplitude  of  the  nth  cycle  and  An+N  is  the  amplitude  of  the  {n+N)th 
cycle.   Similarly  to  a  cantilever,  the  amplitude  of  free  vibration  of  a  whipping  ship  will 
decrease  with  time  due  to  the  presence  of  damping. 

A  second,  and  perhaps  more  important  effect  of  damping  is  on  the  steady  state 
amplitude  attained  when  a  structure  is  excited  by  a  harmonically  oscillating  force.  As 
discussed  in  any  text  on  mechanical  vibrations,  the  maximum  magnitude  of  the  response  of 
any  structure  subjected  to  harmonic  excitation  (or  periodic  excitation  that  can  be 
represented  as  a  Fourier  series  of  harmonic  terms)  is  determined  by  the  relation  of  the 
frequency  of  excitation  to  the  natural  frequencies  of  the  structure,  as  well  as  the  stiffness. 
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Cantilever 


Figure  1-1    Effect  of  damping  (with  time)  on  the  free  vibration 
of  a  cantilever 


damping  and  mass  of  the  structure  (and,  of  course,  the  magnitude  of  the  excitation)    As 
shown  in  figure  1-2  (illustrated  for  a  single  degree  of  freedom  system  subjected  to 
harmonic  excitation),  the  response  is  predominantly  controlled  by  stiffness  when  the 
excitation  frequency  (co)  is  much  less  than  the  natural  frequency  (con)  of  the  structure 
(elastic  forces  dominate  the  equations  of  motion  due  to  the  small  acceleration  and  inertia 
forces).  When  the  excitation  frequency  is  much  greater  than  the  natural  frequency,  the 
response  is  predominantly  controlled  by  mass  (inertia  forces  dominate  the  equations  of 
motion)    When  the  excitation  frequency  is  close  to  the  natural  frequency  (near 
resonance),  however,  the  elastic  and  inertia  forces  approximately  balance,  and  the  major 
demand  on  the  excitation  forces  is  to  overcome  system  damping,  implying  that  the  amount 
of  system  damping  controls  the  dynamic  response    In  this  region,  when  damping  is  low, 
the  response  magnitude  is  high,  and  when  damping  is  high,  the  response  amplitude  is  low 
Thus,  the  response  near  resonance  can  be  very  much  greater  than  the  static  response  (co 
approaching  zero),  particularly  if  the  damping  is  very  low. 

As  discussed  previously,  the  bubble  pulse  can  generate  low  frequency,  quasi- 
periodic  loading  on  the  ship.  For  certain  underwater  explosive  combinations  (explosive 
charge  weight  and  standoff),  the  frequency  of  this  loading  can  be  on  the  order  of  the 
natural  frequencies  of  the  ship's  low  frequency  flexural  modes    Thus,  hull  flexural 
vibrations,  excited  by  the  quasi-periodic  bubble  pulse  loading,  can  approach  a  resonant 
response  in  the  low  frequency  modes.   It  is  clear  then  that  damping  can  have  a  significant 
influence  on  the  magnitude  of  the  whipping  response  of  the  hull  to  the  bubble  pulse 
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Figure  1-2.  Effect  of  damping  (with  frequency)  on  the  steady  state 
vibration  of  a  single  degree  of  freedom  system 
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2.  DAMPING  MECHANISMS  IMPORTANT  IN  SHIP  WHIPPING 

2.1  Introduction. 

In  dynamics,  the  word  damping  generally  refers  to  the  dissipation  of  energy  during 
the  motion  of  a  body  or  structure    This  energy  eventually  is  converted  to  heat  by  friction 
or  generates  waves  in  the  surrounding  medium     When  a  structure,  such  as  a  ship,  vibrates 
or  moves  in  on  ocean  environment,  energy  can  be  dissipated  within  the  structure  itself, 
through  internal  components,  or  externally  into  the  surrounding  fluid     Specific 
mechanisms  of  damping,  which  can  be  expected  to  be  important  in  ship  and  submarine 
whipping,  are  discussed  here  and  then  incorporated  into  a  model  for  submarine  whipping 
in  chapter  3. 

2.2  Hull  damping. 

The  loss  of  vibrational  energy  within  a  complex  structure  can  be  accomplished 
through  several  possible  mechanisms,  depending  on  the  frequency  of  the  vibration  and 
geometry  of  the  structure.  In  general,  low  frequency  energy  (consistent  with  ship 
whipping)  is  lost  within  a  structure  through  friction,  which  in  turn  generates  heat  which  is 
conducted,  convected  and  radiated  away    This  friction  can  occur  among  the  material 
particles  or  crystals  (within  the  structural  elements)  as  they  move  relative  to  one  another 
during  material  deformation    This  damping  mechanism  is  referred  to  as  material  or 
hysteresis  damping.  Friction  can  also  occur  at  the  interface  between  structural  elements 
as  they  slide  relative  to  one  another  (i.e  at  joints).  This  damping  mechanism  is  often 
referred  to  as  structural  or  dry  friction  damping 

As  will  be  discussed  in  greater  detail  subsequently,  it  is  useful  to  define  the  hull 
damping  in  terms  of  equivalent  viscous  damping.  By  doing  this,  the  damping  forces 
associated  with  material  and  structure  can  be  represented  as  being  proportional  to  the 
velocity    The  resulting  equations  of  motion  for  the  vibration  of  the  hull  could  then  be 
written  in  the  form 

Mx  +  Cx  +  Kx  =  f  (2-1) 

where  M  is  the  mass  matrix  of  the  hull,  C  is  the  equivalent  viscous  damping  matrix,  K  is 
the  stiffness  matrix,  x  is  the  displacement  vector,  f  is  the  external  load  vector,  and 
superposed  dots  refer  to  time  derivatives  (x  is  velocity  and  x  is  acceleration).  Although 
material  damping  is  more  representative  of  the  energy  dissipation  within  the  material,  and 
structural  damping  more  representative  of  losses  in  joints,  it  will  be  noted  that  an 
equivalent  viscous  damping  is  assumed  for  most  real  structural  dynamic  analyses.  The 
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precise  nature  of  the  damping  is  less  important  than  modeling  the  correct  energy  loss  per 
cycle 

2.2.1     Material  (hysteresis)  damping. 

When  materials  are  deformed,  energy  is  absorbed  and  dissipated  by  the  material 
due  to  friction  during  internal  reconstruction  of  the  micro  and/or  macro  structure  as  the 
material  deforms.  This  reconstruction  ranges  from  crystal  lattice  to  molecular  scale  effects 
This  process,  referred  to  as  material  damping  or  hysteresis  damping,  is  discussed  in  great 
detail  by  Lazan3  and  Nashif4 

All  real  materials  dissipate  energy  during  cyclic  deformation    Regardless  of  the 
precise  physical  mechanisms  involved,  energy  dissipation  can  be  highly  nonlinear,  and 
therefore  detailed  analysis  can  be  very  difficult  One  approach  toward  quantifying  the 
internal  damping  behavior  of  real  materials  is  through  the  hysteresis  loop,  figure  2-1    The 
loop  plots  strain  vs.  stress  during  cyclic  deformation  of  a  material  volume  and  is  obtained 
from  experimental  measurements.  Typically,  hysteresis  loops  representative  of 
construction  metals  such  as  steel  are  extremely  thin,  deviating  little  from  a  single  line 
(unless  the  metal  is  strained  well  into  the  plastic  range)    Thus,  for  elastic  ship  hull 
whipping  of  a  steel  hull,  the  damping  provided  by  material  hysteresis  would  be  expected  to 
be  small  compared  to  other  damping  mechanisms  discussed  subsequently. 

The  area  inside  the  hysteresis  loop  is  the  specific  damping  energy,  D,  and  is  equal 
to  the  energy  dissipated  per  unit  volume  of  material  per  cycle 

ds 
dt 

where  o  is  stress  and  e  is  strain.  Equation  (2-2)  implies  that  the  rate  of  energy  dissipation 
(dD/dt)  is  proportional  to  the  strain  rate.  This  is  equivalent  to  the  hysteresis  damping  being 
proportional  to  the  strain  rate. 

The  specific  damping  energy  is  very  small  for  most  conventional  structural 
materials,  such  as  steel,  for  typical  elastic  stress-strain  levels.  High-damping  alloys  and 
viscoelastic  materials  have  higher  specific  damping  energies,  and  therefore  are  much  better 
at  dissipating  energy. 


D  =  f<7.«fr=    Ja-    -  U  (2-2) 


3Lazan.  B.J..  Damping  of  Materials  and  Members  in  Structural  Mechanics,  Pergamon  Press.  New  York. 

1968. 

4Nashif.  AD.  Jones.  DIG.  and  Henderson.  J. P..  Vibration  Damping.  John  Wiley  and  Sons.  New  York. 

1985. 
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slightly  nonlinear 

highly   nonlinear 


^     strain 


Figure  2-1 .   Typical  hysteresis  loops. 

As  discussed  in  the  previous  section,  it  is  useful  to  represent  hysteresis  damping  in 
terms  of  an  equivalent  viscous  damping    An  equivalent  viscous  damping  coefficient  for 
hysteresis  damping  can  be  found  by  considering  the  harmonic  motion  of  an  equivalent 
spring-viscous  damper  system  (figure  2-2).  For  this  system,  the  force  F(t)  necessary  to 
cause  a  displacement  x(t)  is  given  by 

F(t)  =  kx(t)  +  cx(t)  (2-3) 

where  k  is  the  spring  constant,  c  is  the  viscous  damping  constant,  and  x  is  the  velocity 
(dx/dt)    For  harmonic  motion  of  frequency  co  and  amplitude  X,  the  displacement  can  be 
written  in  the  form 

x(t)  =  Xsin(ot)  (2-4) 

Inserting  equation  (2-4)  into  equation  (2-3)  gives 

F(t)  =  kXsin(ot)  +  cX(ocos(cot)  (2-5a) 


F(t)  =  kx±ccoVx:  -(Xsincot): 
F(t)  =  kxtccoVx^-x2 


(2-5b) 
(2-5c) 


When  F(t)  is  plotted  against  x(t),  equation  (2-5c)  represents  a  closed  loop  similar  to  a 
hysteresis  loop  (figure  2-3)    The  energy  dissipated  by  the  damper  in  one  cycle  of  motion 
is  found  by  substituting  equation  (2-5a)  into  equation  (2-2) 
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[J     \  -na 

-   dt  =    f  (kXsincot  +  cXcocoscot)(coXcoscot)dt  =7tcocX: 


(2-6) 


where  AW  is  the  energy  loss  in  the  damper  in  one  cycle,  which  is  the  equivalent  to  the 
specific  damping  energy,  D,  for  hysteretic  damping. 

/////////// 


k 


Tx 


(t) 


FCt) 


Figure  2-2.  Equivalent  spring-viscous  damper  system. 

F 


Figure  2-3.  Equivalent  spring- viscous  damper  system 


For  real  hysteresis  damping,  however,  it  has  been  found  experimentally  that  the 
specific  damping  energy,  D,  is  independent  of  frequency,  but  approximately  proportional 
to  the  square  of  the  amplitude  of  the  strain5    Thus,  the  equivalent  viscous  damping 
coefficient,  ceq,  must  be  inversely  proportional  to  the  frequency,  which  can  be  represented 
by 


CeM     = 


CO 


(2-7) 


5Rao.  S ..  Mechanical  Vibrations.  Addison-Welleslev  Publishing  Co..  Reading.  MA.  1990.  p.  102 
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where  h  is  called  the  hysteresis  clamping  constant    As  discussed  subsequently,  the  idea 
that  the  equivalent  viscous  damping  constant  for  material  hysteresis  damping  is  inversely 
proportional  to  frequency  is  made  use  of  in  the  practical  modeling  of  material  damping  in 
real  structures. 

2.2.2     Structural  (dry  friction)  damping. 

In  addition  to  damping  within  the  material  of  a  structure,  damping  also  occurs  at 
the  interface  between  structural  elements  as  they  slide  relative  to  one  another    This  type 
of  damping  is  usually  referred  to  as  structural  damping,  slip  damping,  Coulomb  damping, 
or  simply  dry  friction  damping  (it  amounts  to  Coulomb  or  "dry"  friction  forces)    This 
type  of  damping  is  considered  important  when  structural  elements  are  joined  non-rigidly 
by  riveting,  clamping,  or  bolting,  where  a  moderate  amount  of  clamping  pressure  allows 
for  some  slip  between  elements    As  for  material  damping,  structural  damping  is  discussed 
in  great  detail  by  Lazan6,  and  NashiF. 

Because  of  the  potential  geometric  complexities,  quantitative  evaluation  of 
structural  damping  in  complex  structures  is  impractical  even  when  numerical  techniques 
such  as  finite  element  analysis  are  used    However,  most  structurally  significant  joints 
making  up  a  ship  hull  are  rigidly  welded  rather  than  bolted  or  riveted.   Thus,  for  elastic 
ship  hull  whipping,  it  would  be  expected  that  damping  provided  by  structural  dry  friction 
forces  would  be  very  small  compared  to  other  damping  forces  discussed  subsequently    It 
should  be  noted,  however,  that  for  certain  ship  designs  (particularly  many  submarines  and 
submersibles),  internal  decking  is  often  connected  to  the  main  structure  of  the  hull  through 
various  types  of  sliding  deck  joints  or  suspension  mechanisms    Thus,  in  some  ship 
applications,  it  is  necessary  to  account  for  some  degree  of  dry  friction  damping, 
particularly  when  considering  potentially  large  amplitude  hull  flexural  motions. 

In  general,  Coulomb's  Law  of  Dry  Friction  states  that  when  two  bodies  are  in 
contact,  the  tangential  force  required  to  produce  sliding  is  proportional  to  the  normal 
force  acting  in  the  plane  of  contact.  Once  sliding  has  begun,  the  force  remains  constant. 
The  dry  friction  force  is  given  by 

fD  =  uN  (2-8) 

where  N  is  the  normal  force  and  u  is  the  coefficient  of  friction  which  is  characteristic  for 
the  interface  of  the  two  surfaces.  Thus,  the  structural  dry  friction  damping  force  at  the 
interface  of  two  sliding  surfaces  is  in  a  direction  opposite  to  the  direction  of  the  relative 


6Lazan.  op.  at. 
7Nashif.  op.  cit. 
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velocities  of  the  two  surfaces  but  is  independent  of  the  displacement  or  velocity;  it 
depends  only  on  the  normal  force  between  the  two  surfaces 

The  fact  that  the  dry  friction  forces  remain  constant  as  long  as  sliding  motion  exists 
between  two  surfaces  (and  are  not  dependent  on  displacement  or  velocity)  leads  to  the 
conclusion  that,  in  each  successive  cycle  (or  oscillation)  of  a  structure,  the  amplitude  of 
the  motion  is  reduced  by  an  amount  proportional  only  to  the  friction  forces  resulting  from 
the  sliding  within  the  structure    For  a  simple  one  dimensional  case,  the  reduction  in 
amplitude  with  each  successive  oscillation  will  remain  constant  as  long  as  relative  motion 
exists. 

2.2.3     Modeling  material  and  structural  damping  in  complex  structures. 

If  the  material  and  structural  damping  properties  of  a  structure  could  be 
quantitatively  determined  throughout  the  structure,  a  finite  element  procedure  could  be 
used  to  include  the  distributed  damping  properties  (as  element  damping  matrices)  in  the 
complete  dynamic  solution  of  the  structural  motions.  In  practice,  however,  direct 
evaluation  of  distributed  structural  damping  properties  in  complex  structures  does  not 
lend  itself  well  to  the  numerical  solutions  for  structural  dynamics  (although  distributed 
damping  has  been  included  in  some  simple  finite  element  applications8).  Material  and 
structural  damping  (as  previously  defined)  are  typically  expressed  in  terms  of  damping 
ratios  established  from  experiments,  rather  than  by  means  of  an  explicit  damping  matrix.  In 
fact,  for  many  applications,  there  is  no  need  to  express  the  damping  explicitly  by  means  of 
a  damping  matrix,  but  rather  only  in  terms  of  the  modal  damping  ratios  (Qn)  which  can  be 
applied  through  the  principle  of  modal  superposition.  This  concept  is  often  called  modal 
damping.  An  extension  of  modal  damping  in  which  an  explicit  damping  matrix  is  defined 
through  a  linear  combination  of  the  mass  and  stiffness  matrices  is  called  Rayleigh 
damping. 

Modal  damping  and  modal  superposition.  The  method  of  dynamic  analysis 
using  modal  superposition  is  discussed  in  numerous  texts  (Clough  and  Penzien9,  Smith10, 
for  example)  and  is  discussed  briefly  here  only  for  continuity    Fundamentally,  the  original 
set  of  N  coupled  equations  of  motion  (for  an  N  degree  of  freedom  system)  are 
transformed  to  a  set  of  N  independent  normal  (modal)  coordinate  equations.  The 


8Kaliske.  M  .  Jagusch.  J ..  Gebbeken.  N  .  and  Rothert.  H..  "Damping  models  in  finite  element 
computations".  Structural  Dynamics-Proceedings  of  the  2nd  European  Conference  on  Structural 
Dynamics  (Vol.  I).  Rotterdam.  1993.  pp.  585-591. 

9Clough.  R.W  and  Penzien.  J..  Dynamics  of  Structures  (second  edition).  McGraw-Hill.  New  York.  1993. 
10Smith.  J.W..  I  'ibration  of  Structures  -  Applications  in  Civil  Engineering  Design.  Chapman  and  Hall. 
London.  1988. 
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dynamic  response  can  be  obtained  by  solving  separately  for  the  response  of  each  normal 
(modal)  coordinate  and  then  superposing  these  to  obtain  the  response  in  the  original 
coordinate  system 

The  equations  of  motion  for  a  multiple  degree  of  freedom  dynamic  systems  can  be 
written  in  matrix  form: 

Mx  +  Cx  +  Kx  =  f  (2-9) 

where,  in  general,  matrices  M,  C,  and  K  all  have  non-diagonal  terms  (i.e.  the  equations  of 
motion  are  coupled)    The  modal  equations  of  motion  are  obtained  from  equations  (2-9) 
by  applying  a  modal  transformation, 

x  =  Oq  (2-10) 

where  <X>  is  the  matrix  whose  columns  are  the  mode  shape  vectors 

*  =  [*i4>2  ■••♦„]  (2-H) 

((J)j  is  the  vector  representing  the  shape  of  the  /th  mode),  and  q  is  a  vector  of  normal  or 
modal  coordinates.  Pre-multiplying  each  term  by  the  transform  of  the  nth  mode  shape 
vector,  yields  an  equation  of  motion  for  each  mode  (uncoupled  from  other  modes)  of  the 
form: 

Mnqn+Cnqn+Knqn=Qn  (2-12) 

where: 

Mn  =  4>nrM(J)n  =  modal  mass  coefficient 

Cn  =  4>JC(t)n  =  modal  viscous  damping  coefficient 

Kn  =  (j)„K(J)n  =  modal  stiffness  coefficient 

Qn  =  $TJ  =  modal  force  (2-13) 

Equations  (2-13)  for  modal  mass  and  modal  stiffness  come  about  because  of  the 
orthogonality  property  of  the  normal  mode  shapes,  which  can  be  extended  to  modal 
viscous  damping  if  the  orthogonality  condition  is  also  assumed  to  apply  to  the  damping 
matrix  (as  will  be  discussed  subsequently).  If  equation  (2-12)  is  divided  by  the  modal 
mass  coefficient,  the  modal  equations  of  motion  may  be  written  in  an  alternate  form: 

qn+2Cnconqn+co;qn=^  (2-14) 

where  co  n  is  the  (modal)  natural  frequency,  and  C,n  is  defined  as  the  modal  viscous 
damping  ratio 

^=v^tr  (2-15) 

2©.M. 
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The  (modal)  natural  frequencies  (co  n)  and  mode  shape  vectors  (<j)n)  are  found  by  solving 
the  generalized  eigenproblem  for  the  undamped  free  vibration  of  the  dynamic  system: 

K<t)n=co^M(J)n  (2-16) 

As  discussed  previously,  it  is  often  more  convenient  (and  physically  reasonable)  to 
define  the  damping  of  a  complex  multiple  degree  of  freedom  system  using  a  viscous 
damping  ratio  for  each  mode,  rather  than  to  try  to  explicitly  evaluate  coefficients  of  the 
damping  matrix  C,  because  the  modal  damping  ratios  can  be  determined  experimentally 
(or  estimated  with  adequate  precision  in  many  cases). 

With  the  (uncoupled)  modal  equations  of  motion  defined,  the  modal  displacements 
can  be  found,  and  the  total  displacement  can  be  found  by  summing  the  modal 
contributions  using  equation  (2-10): 

n=l 

Thus,  the  total  response  of  a  system  modeled  using  modal  superposition  can  be  thought  of 
as  a  sum  of  the  responses  oi  N  independent  single  degree  of  freedom  systems. 

Rayleigh  (proportional)  damping.  Under  some  circumstances,  it  is  desirable  or 
practicable  to  use  modal  superposition  to  uncouple  the  equations  of  motion  in  solving  for 
the  dynamic  response  of  a  system.  In  other  circumstances,  it  may  be  desirable  to  use  the 
basic  concept  of  modal  superposition,  but  develop  an  explicit  damping  matrix  (C)  which 
may  be  used  to  solve  the  complete  equations  of  motion  in  the  time  domain.  This  can  be 
done  by  applying  what  has  become  known  as  Rayleigh  damping  or  proportional  damping 
(named  after  Lord  Rayleigh,  who  first  suggested  its  use).  Rayleigh  or  proportional 
damping  is  discussed  in  varying  detail  in  numerous  sources  (Clough  and  Penzien11, 
Bathe12,  for  example). 

Under  Rayleigh  or  proportional  damping,  it  is  assumed  that  the  damping  matrix 
can  be  adequately  represented  by  making  it  a  linear  superposition  of  the  mass  matrix  (mass 
proportional  damping),  the  stiffness  matrix  (stiffness  proportional  damping),  or  both.  By 
applying  a  linear  superposition,  the  orthogonality  condition  which  applies  to  the  mass  and 
stiffness  matrices  (as  discussed  previously)  would  also  apply  to  the  damping  matrix.  Thus, 
the  damping  matrix  could  be  written  in  the  form: 

C  =  ctM  +  pK  (2-18) 


11  Clough  and  Penzien,  op.  cit. 

12Bathe,  K.J.,  Finite  Element  Procedures  in  Engineering  Analysis,  Prentice-Hall,  Englewood  Cliffs,  NJ, 

1982. 
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where  a  and  (3  are  proportionality  constants  The  damping  associated  with  equation  (2-18) 
may  be  recognized  by  evaluating  the  generalized  modal  damping  associated  with  it  by 
combining  it  with  equations  (2-13)  and  (2-15),  and  solving  for  the  modal  damping  ratio 

Cn  =<t>;C<i)n  =4)nT[aM  +  pK](J)n  =aMn  +(3Kn  (2-19) 

Cn  a        Bco  n 

Cn  = —  = +  ^^  (2-20) 

2onMn      2con        2  ' 

n  n  n 

Thus,  the  damping  associated  with  any  particular  mode  can  be  thought  of  as  a  linear 
combination  of  a  factor  which  is  inversely  proportional  to  frequency  (the  mass 
proportional  term),  and  a  factor  which  is  directly  proportional  to  frequency  (the  stiffness 
proportional  term)    An  important  point  to  note  is  that  the  Rayleigh  damping  model  is 
consistent  with  the  discussion  of  material  (hysteresis)  and  structural  damping  presented 
earlier    Thus,  Rayleigh  damping  can  be  thought  of  as  providing  a  good  model  to  account 
for  material  and  structural  damping  (for  elastic  whipping,  as  discussed  in  the  previous 
sections). 

The  relationship  between  damping  ratio  and  frequency  for  Rayleigh  damping  can 
be  illustrated  as  shown  in  figure  2-4.  The  cases  of  purely  mass  proportional  damping 
((3  =  0)  and  purely  stiffness  proportional  damping  (a  =  0)  are  shown  separately  from  the 
combined  or  general  case. 

As  illustrated  in  figure  2-4,  the  two  Rayleigh  damping  proportionality  constants,  a 
and  P,  can  be  evaluated  by  solving  two  simultaneous  equations  if  the  (modal)  damping 
ratios,  Qa  and  £b,  associated  with  two  specific  (modal)  frequencies,  cod  and  co  b,  are 
known  (i.e.  from  experiment).   Substitution  of  the  two  damping  ratios  and  frequencies  into 
equation  (2-20),  combining  the  resulting  equations  into  matrix  form,  and  carrying  out  a 
matrix  inversion  gives  the  solution  of  the  two  damping  proportionality  constants: 


CO  ,03, 

=  2     /    b 


1  /  CO  a       CO 

1/gk     co 


C0h  "IB, 

—  1  /  CO  b      1  /co 


(2-21) 


(2-22) 


cob-co 

It  is  apparent  from  figure  2-4  that  the  use  of  Rayleigh  damping  provides  for  very 
high  damping  ratios  for  modes  of  very  high  frequencies  (co  much  greater  than  co  b).  The 
end  result  is  that  the  responses  of  very  high  frequency  modes  are  effectively  eliminated  by 
the  high  damping  imposed  by  the  Rayleigh  damping.   In  the  case  of  elastic  ship  whipping, 
where  only  the  lower  few  modes  are  of  concern  (as  discussed  previously),  this  limitation 
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can  be  considered  acceptable  as  long  as  the  higher  frequency  used  for  evaluation  of  the 
Rayleigh  proportionality  constants  (o  h)  is  set  among  the  higher  modes  expected  to 
contribute  significantly  to  the  whipping  response. 

Damping   ratio 


(Jo 


Frequency 


Figure  2-4.  Relationship  between  damping  ratio  and  frequency  for  Rayleigh  damping. 

2.3        Internal  damping. 

Within  any  complex  structure,  such  as  a  ship,  energy  can  be  lost  (or  simply 
transferred)  to  "separate"  dynamic  systems  which  are  in  some  form  "attached"  to  the  main 
structure    For  ships,  two  broad  categories  could  be  considered  under  this  category. 
Discrete  internal  masses  or  structures  such  as  machinery  and  equipment  connected 
indirectly  to  the  hull  through  decking,  machinery  foundations,  and  shock  or  sound  mounts, 
can  act  as  dynamic  subsystems,  redistributing  and  absorbing  dynamic  energy  transferred 
from  the  hull.   Additionally,  onboard  liquids  which  have  free  surface  can  act  as  dynamic 
subsystems,  redistributing  and  absorbing  dynamic  energy  by  sloshing.  Because  of  the 
large  number  of  discrete  internal  structures  and  liquid  storage  tanks  found  onboard  most 
ships,  it  is  important  to  consider  the  mechanisms  through  which  energy  can  be  transferred 
and  absorbed  into  these  internal  damping  devices. 

2.3.1     Damping  through  discrete  internal  structures. 

Heavy  machinery,  electronic  equipment,  and  other  discrete  masses  onboard  a  ship 
can  be  considered  to  have  a  significant  potential  for  dissipating  hull  whipping  energy. 
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Figure  2-5  illustrates  notional  discrete  masses,  connected  to  the  hull  via  mounts,  internal 
decking,  foundations,  etc.  As  will  be  discussed  subsequently,  the  dynamic  response  of 
more  complex  internal  structures  could  be  modeled  as  easily  as  individual  masses  by 
applying  the  principle  of  modal  superposition  (under  the  assumption  of  linear  response). 
For  purposes  of  this  analysis,  the  term  "discrete  internal  structure"  will  be  used  to  refer  to 
any  mass  or  discrete  structure  that  may  be  modeled  as  a  discrete  internal  dynamic  system 
in  terms  of  transfer  and  absorption  of  hull  whipping  energy    While  it  is  certainly  not 
desired  that  heavy  machinery  or  electronic  equipment  absorb  large  amounts  of  energy,  it 
must  be  assumed  that  their  ultimate  connection  to  the  hull  provides  the  potential  for  the 
transfer  of  energy  from  the  whipping  hull 


Figure  2-5.  Discrete  masses  connected  to  the  hull  via 
mounts,  internal  decking,  foundations,  etc. 

Whether  or  not  discrete  internal  structures  connected  to  a  hull  provide  damping  to 
hull  whipping  depends  upon  the  sizes  of  the  masses,  the  dynamic  properties  of  the 
mechanisms  through  which  they  are  connected  to  the  hull  (decking,  foundations,  mounts, 
etc.),  the  dynamic  properties  of  the  hull,  and  the  frequency  of  hull  whipping  (frequency  of 
excitation).  In  section  1.2,  it  was  stated  that  the  controlling  mechanism  in  the  response  of 
a  system  to  excitation  (stiffness,  damping,  or  mass)  is  determined  by  the  relation  of  the 
excitation  frequency  (or  frequencies)  to  the  natural  frequency  (or  frequencies)  of  the 
system.  The  maximum  amount  of  energy  is  transferred  when  the  frequency  of  excitation  is 
equal  to  a  natural  frequency  (at  resonance)    In  the  case  of  discrete  internal  structures 
connected  to  a  whipping  hull,  it  is  evident  that  maximum  energy  would  be  transferred 
when  the  frequency  of  hull  whipping  is  equal  to  a  natural  frequency  of  the  system 
composed  of  the  hull,  the  discrete  internal  structure,  and  the  connection  mechanisms.  An 
understanding  of  the  dynamic  effects  of  discrete  internal  structures  connected  to  a 
whipping  ship  can  be  made  by  first  considering  the  effects  of  a  single  mass  connected  to  a 
(rigid)  vibrating  structure,  where  the  structure  and  mass  are  free  to  vibrate  in  only  one 
direction  (a  two  degree  of  freedom  system),  and  then,  using  the  principle  of  modal 
superposition,  extending  this  to  more  complex  arrangements  of  multiple  masses  connected 
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to  a  vibrating  structure,  such  as  a  whipping  ship  hull  (i.e.  multiple  degree  of  freedom 
systems) 

The  two  degree  of  freedom  system  (dynamic  vibration  absorber)    A  single 
mass  connected  to  a  much  larger  vibrating  structure  can  be  used  to  dissipate  or  absorb 
vibrational  energy  of  the  larger  structure    In  machinery  and  structural  dynamics,  this 
concept  has  gained  practical  significance  in  the  application  of  the  dynamic  vibration 
absorber    The  concept  of  the  dynamic  vibration  absorber  is  presented  by  Den  Hartog13, 
and  their  design  and  application  is  discussed  in  detail  in  numerous  other  sources 
(Korenev14,  Hunt15,  for  example)    A  small  auxiliary  mass  is  connected  to  a  harmonically 
vibrating  structure  by  a  spring  of  specified  stiffness  (the  location  of  the  mass  is  selected  to 
have  the  maximum  effect  on  the  significant  resonant  mode  of  the  vibrating  structure)    In 
the  case  where  damping  is  provided  by  the  connecting  mechanism  (e.g.  by  a  dashpot),  the 
auxiliary  mass  is  called  a  damped  dynamic  vibration  absorber  or  tuned  mass  damper. 
Treating  the  main  structure  and  vibration  absorber  as  a  two  degree  of  freedom  dynamic 
system  (figure  2-6),  Den  Hartog  showed  that  the  amplitude  of  harmonic  vibration  of  the 
main  structure  could  be  reduced  to  zero  at  a  specified  excitation  frequency  (or  reduced 
over  a  range  of  frequencies)    This  could  be  accomplished  by  "tuning"  the  properties  of 
the  vibration  absorber  (adjusting  the  mass,  stiffness  of  the  spring,  and  damping  in  the 
dashpot).  The  downside  in  the  application  of  dynamic  vibration  absorbers  appears  in  the 
potentially  large  amplitude  response  (displacement  and  acceleration)  of  the  absorber 
masses    In  the  case  of  heavy  machinery  and  electronic  equipment  onboard  a  whipping 
ship,  this  could  translate  to  high  displacements  and  accelerations  on  the  delicate 
equipment,  and  large  forces  on  the  mounts  and  decking. 

For  the  simple  dynamic  vibration  absorber  attached  to  a  large  vibrating  structure, 
the  equations  of  motion  for  this  simple  two  degree  of  freedom  system  can  be  derived  by 
considering  dynamic  equilibrium  of  each  of  the  masses  separately.  Figure  2-7  illustrates 
each  of  the  masses  and  the  dynamic  forces  acting  upon  them.  From  equilibrium,  the 
equations  of  motion  may  be  written: 

MX  +  KX  =  Fext  +  c(x  -X)  +  k(x  -  X)    (for  the  main  mass)  (2-23) 

mx  +  c(x-  X)  +  k(x-  X)  =  0    (for  the  absorber  mass)  (2-24) 

where  capital  letters  refer  to  the  main  structure  and  lower  case  letters  refer  to  the 
absorber,  as  shown  in  figure  2-6.  These  equations  of  motion  could  be  written  in  an 
alternate  form  by  defining  the  force  exerted  on  the  main  structure  by  the  absorber. 


13Den  Hartog.  J.P ..  Mechanical  Vibrations,  (4th  edition).  McGraw-Hill.  New  York,  1956. 
14Korenev  and  Reznikov,  Dynamic  Vibration  Absorbers.  John  Wiley  and  Sons.  New  York.  1993. 
15Hunt.  J.B..  Dynamic  Vibration  Absorbers.  Mechanical  Engineering  Publications.  London.  1979. 
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absorber 


k(x-X)  +  c(x-X) 


(2-25) 


which,  of  course  is  the  negative  of  the  force  exerted  on  the  absorber  by  the  main  structure 
Thus,  the  equations  of  motion  could  be  written 


MX  +  CX  +  KX  =  F  .  +  F 


abs 


(for  the  main  mass) 


absorber 


=  c(x  -  X)  +  k(x  -  X)  =  -mx    (for  the  absorber  mass) 


(2-26) 

(2-27) 


Equations  (2-26)  and  (2-27)  form  a  set  of  2  simultaneous  differential  equations  (coupled) 
which  may  be  solved  under  various  conditions  for  displacements,  velocities  and 
accelerations. 
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Figure  2-6.  Dynamic  vibration  absorber 
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Figure  2-7    Dynamic  forces  on  each  mass  of  a  dynamic  vibration  absorber 


Multiple  degree  of  freedom  systems.  For  more  complex  arrangements  of 
internal  structure,  where  many  degrees  of  freedom  would  be  required  to  properly  model 
the  dynamic  response,  the  displacements,  accelerations,  and  velocities  would  be  vectors, 
and  a  larger  number  of  simultaneous  differential  equations  would  be  required.  For 
multiple  degree  of  freedom  main  and  internal  structures,  the  equations  of  motion  could  be 
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written  just  as  equations  (2-23)  and  (2-24),  but  in  matrix  form.  Thus,  there  would  be  one 
simultaneous  equation  for  each  degree  of  freedom. 

An  alternative  approach  for  multiple  degree  of  freedom  systems  (particularly 
where  the  internal  structure  is  complex  and  a  large  number  of  degrees  of  freedom  would 
be  required  to  adequately  represent  the  dynamics),  would  be  to  use  the  principle  of  modal 
superposition  under  the  assumption  of  linear  response  (as  discussed  in  section  2  2  3)  to 
approximate  the  response  of  the  internal  structure    By  including  only  those  resonant 
response  modes  of  the  main  structure  and  the  internal  structure  which  would  contribute 
significantly  to  the  overall  response,  the  complexity  of  approximating  the  response  of  the 
main  structure  is  greatly  reduced    Figure  2-8  depicts  a  notional  complex  internal 
arrangement  of  discrete  masses  with  many  degrees  of  freedom. 


discrete   Internal   structure 


X(t) 


Figure  2-8.  Complex  internal  arrangement  of  discrete  masses  with 
many  degrees  of  freedom 

Using  modal  superposition,  an  equivalent  mechanical  model  of  the  dynamic 
behavior  of  the  discrete  internal  structure  may  be  made  as  shown  in  figure  2-9.    Each 
mode  of  resonant  response  of  the  discrete  internal  structure  is  represented  by  an 
equivalent  modal  mass,  modal  damping  and  modal  stiffness  (mass  mn  represents  that 
portion  of  the  discrete  internal  structure  responding  as  a  rigid  body  with  the  main 
structure).  A  variation  of  the  dynamic  equilibrium  equation  (2-24)  for  the  dynamic 
vibration  absorber  can  be  written  using  the  principle  of  modal  superposition    The  equation 
of  motion  for  the  /rth  mode  of  the  internal  structure  can  be  written: 

mnxn+cn(xn-X)  +  kn(xn-X)  =  0  (2-28) 

where  xn  is  the  nth  modal  displacement  of  the  internal  structure  (equivalent  to  qn  in 
section  2.2.3)  and  X  is  the  displacement  of  the  main  structure  (hull). 

The  equations  of  motion  for  a  system  made  up  of  an  internal  structure  attached  to 
a  large  vibrating  structure  (the  hull)  can  be  derived  by  considering  dynamic  equilibrium  of 
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the  main  structure  and  the  superposed  modal  masses  of  the  internal  structure.  The  forces 
applied  to  the  main  structure  by  the  internal  structure  can  be  obtained  by  modal 
superposition    Figure  2-10  illustrates  the  dynamic  equilibrium  for  the  main  structure. 
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Figure  2-9.  Equivalent  mechanical  model  of  the  dynamic  behavior  of 
a  discrete  internal  structure  (using  modal  superposition) 


Fint 


Figure  2-10.  Dynamic  equilibrium  for  the  main  structure. 


From  dynamic  equilibrium,  the  equation  of  motion  for  the  main  structure  can  be 


written: 


[MX  +  CX  +  KX]  =  Fext+F 


ext  int 


(2-29) 
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where  Fext  is  an  external  force  applied  to  the  main  structure,  and  Fint  is  defined  as  the  force 
exerted  on  the  main  structure  by  the  internal  structure    Fmt  can  be  determined  by 
considering  dynamic  equilibrium  of  the  internal  structure  using  modal  superposition,  as 
shown  in  figure  2-11.   By  dynamic  equilibrium  of  the  internal  structure. 


F*  =  -im„x„^[cJx„-X)  +  kJx„-X) 


(2-30) 


Ek(x-X)  !     £c(*-X) 

Figure  2-11    Dynamic  equilibrium  of  the  internal  structure  using  modal  superposition. 


It  should  be  noted  that  the  internal  structure  rigid  body  case  (n  =  0)  can  be  taken  out  of 
equation  (2-30)  by  noting: 

X[cn(xn-X)  +  kJxn-X)]  =  -Xmnxn=-m0X-Xmnxn=-m0X  +  X[cJxn-X)  +  kn(xn-X)] 

n=0  n=0  n=l  n=l 

(2-31) 
With  this  simplification,  the  force  exerted  on  the  main  structure  by  the  internal  structure 
can  be  written: 


Fte=-m0X  +  2[cn(xB-X)  +  kB(xB-X) 

n  =  l 


(2-32) 


Equations  (2-29),  (2-28),  and  (2-32)  form  a  set  of  N+I  simultaneous  equations  (coupled) 
which  can  be  solved  under  various  conditions  for  displacements,  velocities  and 
accelerations. 

It  should  be  noted  that,  for  cases  where  more  than  one  degree  of  freedom  is 
important  for  the  main  structure  (e.g.  the  main  structure  moves  horizontally  as  well  as 
vertically),  the  above  method  could  be  extended  considering  the  modal  response  of  the 
internal  structure  in  each  of  the  degrees  of  freedom  of  the  main  structure    Thus,  the  main 
structure  displacement  would  be  represented  by  a  vector  (X),  the  internal  structure  modal 
displacement  would  be  represented  by  a  vector  (xn),  and  the  forces  (Fm  and  FexI )  would 
be  represented  by  vectors 
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An  application  of  damping  through  discrete  internal  structures  will  be  included  in 
the  development  of  a  model  for  submarine  whipping  in  chapter  3 
2.3.2     Damping  through  sloshing  liquids. 

Liquid  sloshing  has  been  the  focus  of  considerable  attention  in  the  space  industry, 
due  primarily  to  concerns  about  fuel-rocket  interaction  forces,  and  recognition  that  liquid 
sloshing  resonance  can  result  in  significant  dissipation  of  energy  and  can,  therefore,  be 
capitalized  upon  to  suppress  vibrations  in  large  structures16    Similarly,  methods  to 
account  for  the  effects  of  sloshing  liquids  on  the  natural  frequencies  and  damping  of 
offshore  platforms  have  been  developed  within  the  marine  industry17    The  occurrence  of 
resonance,  effects  of  nonlinearities,  and  contribution  of  devices  such  as  baffles  have  been 
the  object  of  numerous  other  studies    It  can  be  projected  that  liquids  in  tanks  onboard  a 
whipping  ship  could  undergo  sloshing  resonance,  resulting  in  the  damping  of  whipping 
energy    Fundamentally,  any  onboard  liquid  with  free  surface  has  the  potential  for 
absorbing  whipping  energy  through  liquid  sloshing  (liquid  without  free  surface  acts  only  as 
a  rigid  mass).  Figure  2-12  illustrates  a  notional  liquid  storage  tank,  restrained  to  move 
with  the  whipping  hull,  and  the  resulting  sloshing  (or  standing  waves)  generated  due  to  the 
existence  of  the  free  surface  of  the  liquid  in  the  tank. 


Figure  2-12.   Tank  with  sloshing  liquid. 

The  liquid  in  a  moving  tank  can  be  expected  to  respond  in  a  variety  of  ways 
dependent  upon  many  variables  (tank  geometry,  directions  of  tank  motion,  amplitudes  and 
frequencies  of  tank  motion,  properties  of  the  liquid,  etc  )    Translational  or  pitching 


16Welt.  F  ,  and  Modi.  V.J.,  "Vibration  damping  through  liquid  sloshing".  Diagnostics,  Vehicle  Dynamics 
and  Special  Topics:  12th  Biennial  Conference  on  Mechanical  Vibrations  and  Noise.  ASME.  1989. 
17Vandiver.  JK.,  and  Mitome.  S..  "Effect  of  liquid  storage  tanks  on  the  dynamic  response  of  offshore 
platforms".  Dynamic  Analysis  of  Offshore  Structures:  Recent  Developments.  CML  Publications. 
Southampton.  1982. 
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motions  of  a  tank  lead  primarily  to  lateral  sloshing  (usually  antisymmetric  liquid  mode 
shapes)    Tank  motions  normal  to  the  equilibrium  free  surface  lead  primarily  to  vertical 
sloshing  (usually  symmetric  liquid  mode  shapes)    Under  conditions  of  abrupt  or  sudden 
acceleration,  liquid  impact  with  tank  overheads  and  opposing  bulkheads  may  provide 
significant  loading  on  tank  boundaries    The  high  accelerations  associated  with  a  whipping 
ship  may  require  consideration  of  any  or  all  of  these  to  account  for  the  transfer  of 
whipping  energy. 

Analysis  of  liquid  sloshing  can  be  carried  out  by  numerous  means    Researchers  in 
the  field  of  liquid  sloshing  have  used  several  principle  methods  to  determine  the  important 
dynamic  characteristics  desired  for  various  applications18    In  many  applications,  potential 
flow  models,  using  linearized  and  nonlinear  free  surface  conditions,  have  been  developed 
and  used  to  approximate  the  fluid  motion  in  various  container  shapes  (Bauer'9,  Hutton20, 
Welt  and  Modi21,  for  example).  In  order  to  properly  account  for  energy  transfer  and/or 
dissipation,  incorporation  of  liquid  impact  loading  and  viscous  effects  has  been 
accomplished  using  equivalent  mechanical  models  in  conjunction  with  potential  flow 
models  (Dalzell22,  Vandiver  and  Mitome23,  Bauer24,  for  example) 

Potential  flow  model.  Fundamentally,  lateral  and  vertical  sloshing  can  be  thought 
of  as  the  motion  of  standing  waves  within  a  moving  rigid  container  or  tank.   Under  the 
assumptions  of  incompressible  and  irrotational  flow,  the  motion  of  the  liquid  in  a  tank  can 
be  approximated  well  by  solving  the  potential  flow  problem  for  the  wave  motion  subject  to 
appropriate  boundary  conditions  (all  variables  are  functions  of  time).  The  velocity 
potential  function  <t>  (defined  by  v  =  V  O  where  v  is  the  fluid  velocity  vector  and  V  is  the 
vector  differential  operator)  represents  solution  of  the  differential  equation  (Laplace 
equation) 

V:O  =  0    (within  the  fluid  domain)  (2-33) 

subject  to  the  following  boundary  conditions: 


18Abramson.  H.N.  (ed).  The  Dynamic  Behcnnor  of  Liquids  in  Moving  Containers,  NASA  SP-106.  1966. 

19Bauer,  H.F.,  "Theory  of  the  Fluid  Oscillations  in  a  Circular  Ring  Tank  Partially  Filled  with  Liquid", 

NASA  Technical  Note  D-5 57,  1960. 

20Hutton.  RE..  "An  Investigation  of  Resonant.  Nonlinear.  Non-Planar  Free  Surface  Oscillations  of  a 

Fluid",  NASA  Technical  Note  D- 1870.  1963. 

2 'Welt  and  Modi.  op.  cit. 

22Dalzell.  J.F  .  "Liquid  Impact  on  Tank  Bulkheads".  The  Dynamic  Behavior  of  Liquids  in  Moving 

Containers,  NASA  SP-106.  H.N.  Abramson  (ed).  1966.  pp.  353-372. 

23Vandiver  and  Mitome.  op.  cit. 

24Bauer.  H.F..  "Nonlinear  Mechanical  Model  for  the  Description  of  Propellant  Sloshing".  AIAA  Journal. 

Vol.  4.  No.  9.  1966.  pp.  1662-1668. 
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=  Vn    (on  the  wall  of  the  tank)  (2-34a) 

dn 

c:0        rO  r<b      1 

— -  +  g +  2VOV +  -VOV(VOVO)  =  0 

ct~  cy  dt      2 

(at  y  =  r|  (on  the  free  surface))  (2-34b) 

where  n  is  the  normal  vector  of  the  tank  wall,  Vn  is  the  velocity  of  the  tank  wall  normal  to 
its  surface,  g  is  the  acceleration  due  to  gravity,  and  r|  is  the  free  surface  elevation    The 
free  surface  boundary  condition,  equation  (2-34b),  represents  the  complete  nonlinear  free 
surface  boundary  condition25.  The  linearized  free  surface  boundary  condition  (linearized 
about  the  mean  free  surface)  would  be  written: 

^^  +  g— =  0    (aty  =  0)  (2-34c) 

ct"         cy 

The  fluid  velocity  vector  (v)  can  be  obtained  from  (by  definition  of  the  velocity  potential): 

v  =  VO  (2-35) 

The  hydrodynamic  pressure  (Ph)  exerted  on  the  tank  wall  by  the  sloshing  liquid  can  be 

obtained  from  the  unsteady  Bernoulli  equations: 


Ph=-P 


cO      1 , 

+  -  VO 

at      2 


(2-36) 


where  p  is  the  liquid  density.   Finally,  the  hydrodynamic  force  exerted  on  the  tank  wall  by 
the  sloshing  liquid  can  be  obtained  by  integration  of  the  hydrodynamic  pressure  over  the 
area  of  the  tank  wall: 

Fh  =JJ  (Phn)dS  (2-37) 

s 

where  F  is  the  hydrodynamic  force  vector  (on  the  tank  wall)  and  S  is  the  surface  of  the 
tank  wall  (below  the  free  surface).  Thus,  for  a  specified  tank  geometry  and  motion,  the 
hydrodynamic  forces  exerted  by  the  sloshing  liquid  on  the  tank  can  be  approximated  using 
a  potential  flow  model 

Equivalent  mechanical  model.  As  mentioned  previously,  in  order  to  properly 
account  for  energy  dissipation,  viscous  effects  at  the  tank  wall  must  be  included    While  it 
is  difficult  to  include  viscous  effects  in  the  potential  flow  solution,  viscous  damping  may 
be  included  in  an  equivalent  mechanical  model  for  the  sloshing.  In  an  equivalent 


25Ne\vman.  J.N..  Marine  Hydrodynamics.  The  MIT  Press.  Cambridge.  MA.  1986.  p.  247. 
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mechanical  model,  each  of  the  sloshing  modes  (or  modes  of  standing  wave  motion)  is 
represented  by  an  equivalent  spring-mass-damper  as  shown  in  figure  2-13    In  typical 
practice,  linear  viscous  damping  is  usually  assumed  for  the  (modal)  dampers  in  the 
mechanical  model26,  however,  nonlinear  mechanical  modeling  techniques  do  exist27 
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Figure  2-13.  Equivalent  mechanical  model  of  a  tank  with  sloshing  liquid 
using  modal  superposition. 


In  an  equivalent  mechanical  model,  the  size  of  the  equivalent  masses  for  each  of 
the  resonant  sloshing  modes  (m,  ,m: ,...)  may  be  determined  by  requiring  that  the 
(hydrodynamic)  force  exerted  on  the  tank  wall  by  the  equivalent  mass-spring  system,  for 
the  case  of  zero  damping,  be  equal  to  the  force  calculated  from  the  potential  flow  theory 
for  each  resonant  sloshing  mode  (m0  is  the  portion  of  the  liquid  in  the  tank  which  acts  as  a 
rigid  body). 


mn*n  =>mn  = 


x. 


(2-38) 


26Silverman,  S..  and  Abramson.  H.N..  "Damping  of  Liquid  Motions  and  Lateral  Sloshing",  The  Dynamic 
Behavior  of  Liquids  in  Moving  Containers,  NASA  SP-106,  H.N.  Abramson  (ed.).  1966.  pp.  105-144. 
27Bauer.  op.  at. 
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where  fnh  is  the  hydrodynamic  force  exerted  on  the  tank  wall  by  the  mh  sloshing  mode 
The  equivalent  modal  stiffnesses  (k,  ,k: , ...)  may  be  determined  from  the  natural 
frequencies  (co  n)  calculated  from  the  potential  flow  theory  for  each  resonant  sloshing 
mode: 


®n  =J—  =>kn=G);mn  (2-39) 


Damping  coefficients  for  each  resonant  sloshing  mode  can  be  determined 
experimentally  for  a  given  tank  configuration  via  a  logarithmic  decrement  method  (section 
1.2).  The  modal  damping  ratio  for  the  mh  mode  (Qn)  can  be  determined  from  the  log 
decrement  (5)  measured  for  each  mode  by 

C,n  =    .  =  «  —  (for  small  damping)  (2-40) 

V(2tu):+5:       2ti 

where  the  log  decrement  is  defined  by  equation  (1-1).  The  modal  damping  ratio  is  the 
ratio  of  the  modal  damping  coefficient  to  the  critical  modal  damping  coefficient  (c  /  cc), 
where  the  critical  modal  damping  coefficient  represents  that  value  of  damping  where  the 
motion  would  not  oscillate,  but  would  be  non-periodic.  The  modal  damping  coefficients 
(c,  ,c: , ...)  can  be  determined  from 

c    =2Cmncon  (2-41) 

n  ~  n        n       n  v  / 

Referring  to  figure  2-13,  the  equation  of  motion  for  the  mh  equivalent  mass  can  be 
written  in  terms  of  the  displacement,  velocity,  and  the  acceleration  of  the  equivalent  mass 
as  follows: 

mnxn+cn(xn-X)  +  kn(xn-X)  =  0  (2-42) 

where  X  is  the  displacement  of  the  tank,  xn  is  the  displacement  of  the  mh  sloshing  mode 
(equivalent  mass),  and  velocities  and  accelerations  are  represented  by  superposed  dots. 
The  total  hydrodynamic  force  on  the  wall  of  the  tank  can  be  found  by  summing  the  forces 
contributed  by  each  mode: 

Fh=-f>nxn=-m0X-f>nxn  (2-43) 

n=0  n=l 

Similar  to  the  discrete  internal  structure,  each  of  the  resonant  modes  of  liquid  sloshing 
can  be  applied  to  the  main  structure  (the  hull). 

As  an  example,  consider  the  lateral  motion  of  a  rigid  structure  with  an  integral 
laterally  sloshing  tank.  The  equivalent  mechanical  model  could  be  as  shown  in  figure 
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2-14    The  equation  of  motion  for  each  of  the  equivalent  liquid  masses  (modal  masses)  is 
given  by  equation  (2-42),  and  the  total  hydrodynamic  force  on  the  wall  of  the  tank  due  to 
the  liquid  sloshing  is  given  by  equation  (2-43)    The  equation  of  motion  for  the  main 
structure  can  then  be  found  by  adding  the  total  hydrodynamic  force  into  the  equation  of 
motion  for  the  structure  without  the  liquid: 

[MX  +  CX  +  KX]  =  Fext  +  Fh  (2-44) 

where  Fext  is  an  external  force  applied  to  the  main  structure.  The  hydrodynamic  force  (Fh) 
may  be  written  (following  the  same  principles  of  modal  superposition  used  for  the  internal 
structure): 


Fh=-Zrnnxn=-mnX  +  2[cn(^-X)  +  kn(xn-X)] 


(2-45) 


n    ■< 


It  can  also  be  noted  that  m„,  the  portion  of  the  liquid  which  acts  as  a  rigid  body,  can  be 
found  from: 


Ml    ^Z"1"    ^m"    =MI    ~Zmn 


(2-46) 


n=l 


where  M,  is  the  total  mass  of  liquid  in  the  tank. 


Figure  2-14.  Equivalent  mechanical  model  for  the  lateral  motion  of  a  rigid 
structure  with  an  integral  laterally  sloshing  tank. 
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In  terms  of  real  applications,  such  as  in  the  whipping  of  a  ship  hull,  equation  (2-45) 
may  be  truncated  because  the  equivalent  masses  mn  and  modal  displacements  xn   grow 
smaller  with  n  (the  rate  at  which  they  grow  smaller  depends  upon  the  geometry  of  the 
tank).  AH  terms  should  be  included  for  which  the  natural  frequencies  of  the  tank  (co  n )  are 
less  than  or  equal  to  the  natural  flexural  frequency  of  the  overall  structure  (less  than  the 
fundamental  flexural/whipping  natural  frequency).   If  equation  (2-45)  is  truncated  to 
include  TV  modes,  then  equations  (2-44),  (2-45),  and  (2-46)  form  a  set  of  TV-  / 
simultaneous  equations  of  motion  for  the  total  dynamic  system  consisting  of  the  main 
structure  and  the  sloshing  liquid. 

It  should  be  noted  that,  for  cases  where  more  than  one  degree  of  freedom  is 
important  for  the  main  structure  (e.g.  the  main  structure  moves  horizontally  as  well  as 
vertically),  the  above  method  could  be  extended  considering  the  modal  response  of  the 
sloshing  liquid  in  each  of  the  degrees  of  freedom  of  the  main  structure    Thus,  the  main 
structure  displacement  would  be  represented  by  a  vector  (X),  the  sloshing  liquid  modal 
displacement  would  be  represented  by  a  vector  (xn),  and  the  forces  (Fh  and  F,xt )  would  be 
represented  by  vectors. 

An  application  of  liquid  sloshing  will  be  included  in  the  development  of  a  model  for 
submarine  whipping  in  chapter  3. 

2.4        External  (fluid)  damping. 

The  loss  of  vibrational  energy  of  a  whipping  ship  into  the  surrounding  fluid 
medium  can  be  accomplished  through  several  general  mechanisms.  If  the  ship  is 
sufficiently  close  to  the  free  surface,  or  its  motion  is  of  sufficient  velocity  and/or 
frequency,  surface  and  pressure  waves  of  non-negligible  energy  can  be  generated.  This 
form  of  fluid  damping  is  often  referred  to  as  wave  radiation  damping.  Additionally, 
whipping  energy  can  be  lost  through  mechanisms  associated  with  the  viscosity  of  the 
surrounding  fluid.  These  mechanisms  can  be  lumped  under  the  general  heading  of  viscous 
fluid  damping,  but  include  specific  mechanisms  associated  with  time-dependent  viscous 
boundary  layer  shear  forces  ("skin  friction"  drag)  as  well  as  boundary  layer  separation 
("form"  drag). 

In  general,  the  dominant  hydrodynamic  forces  on  a  structure  are  determined 
primarily  by  the  characteristics  of  the  fluid  flow  around  the  structure    For  a  whipping  ship, 
the  relative  fluid  flow  around  the  hull  would  be  generally  oscillatory  in  nature.  For 
oscillatory  flow,  the  following  non-dimensional  parameters  are  often  used  to  characterize 
the  flow: 
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Peak  Reynolds'  number,  Rn  =  — - — 


U  T 

Keulegan-Carpenter  number,  KC  =  — — 

Rn       D; 


Viscous-frequency  parameter,  J3  = 


KC      vT 

where  U  is  the  characteristic  relative  free  stream  velocity  (U,,  is  the  amplitude  of  the 

characteristic  relative  free  stream  velocity  or  maximum  relative  velocity),  D  is  the 
characteristic  dimension  of  the  structure  (i.e.  the  diameter),  v  is  the  kinematic  viscosity  of 
the  fluid,  and  T  is  the  characteristic  period  of  oscillation  of  the  relative  fluid  flow 
Additionally,  the  flow  can  be  influenced  significantly  by  the  geometry  of  the  structure,  free 
surface  effects,  sea  floor  effects,  roughness  on  the  surface  of  the  structure,  and  the 
orientation  of  the  structure  relative  to  the  flow    It  may  also  be  shown  that  for 
harmonically  oscillating  flow  around  a  fixed  circular  cylinder,  KC  =  2tiA  /  D,  where  A  is 
the  amplitude  of  the  oscillation  of  the  fluid  particle  displacement.  Thus,  KC  may  be 
thought  of  as  representing  an  amplitude  of  motion  between  a  structure  and  the  fluid 
relative  the  characteristic  dimension  of  the  body. 

These  parameters,  in  addition  to  the  other  influencing  factors  listed  above, 
determine  the  characteristics  of  the  flow,  and  thus  the  major  components  of  the 
hydrodynamic  force  on  the  structure.  As  will  be  discussed  in  greater  detail  subsequently, 
these  characteristic  parameters  (i.e.  the  characteristics  of  the  flow)  effect  the  occurrence  of 
boundary  layer  separation,  and  thus  a  major  component  of  viscous  fluid  forces  (as  will  be 
discussed  subsequently,  boundary  layer  separation  forces  tend  to  dominate  over  boundary 
layer  shear  forces  for  bluff  bodies)    For  fully  submerged  bodies,  if  KC  is  large  (i.e.  U0  is 
large  and  D  is  small),  then  boundary  layer  separation  is  likely  to  occur,  and  the  major 
hydrodynamic  forces  are  likely  to  be  associated  with  viscous  drag  effects  (the  so-called 
"drag-dominated  regime")      If  both  Rn  and  KC  are  large  (i.e.  U0  is  large  but  D  is  not 
small),  then  the  added  forces  associated  with  the  inertia  of  the  surrounding  fluid  (called 
fluid  inertia  forces)  must  also  be  considered  in  addition  to  the  drag  forces  (the  so-called 
"Morison  regime")    If  KC  is  small  (Rn  may  still  be  moderately  large),  then  boundary 
layer  separation  is  not  likely  to  occur  and  the  only  viscous  forces  would  be  due  to 
boundary  layer  shear    If,  in  addition  to  KC  being  small,  the  structure  is  on  or  near  a  free 
surface  then  radiation  of  surface  waves  tends  to  dominate  the  hydrodynamic  damping  (the 
so-called  "diffraction/radiation"  regime).  It  is  important  to  note  that,  for  complex 
structures  such  as  ships,  different  portions  of  the  structure  may  be  subject  to  different 
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types  of  hydrodynamic  loading  (depending  upon  the  local  characteristics  of  the  structure 
and  the  relative  fluid  flow). 


2.4.1     Wave  radiation  damping. 

Fundamentally,  waves  are  radiated  whenever  a  structure  moves  in  a  fluid  medium 
As  mentioned  above,  when  KC  is  small  (i.e.  when  D  is  large),  and  the  structure  is  on  or 
near  a  free  surface,  then  forces  associated  with  wave  radiation  are  generally  important    In 
typical  practice,  panel  or  boundary  element  methods  are  used  to  solve  the  boundary-value 
problem  in  analyzing  hydrodynamic  forces  due  to  incident,  diffracted  and  radiated 
pressures    Boundary  element  methods  are  based  upon  potential  flow  theory,  which 
neglects  the  effects  of  viscosity.  However,  there  has  been  some  effort  in  recent  years  to 
include  some  of  the  effects  of  viscosity  within  the  boundary  element  method  (as  will  be 
discussed  subsequently) 

Potential  flow  theory  and  its  application  to  the  hydrodynamics  of  structures  is 
discussed  in  great  detail  in  numerous  sources  (Newman28,  Chakrabarti29,  Faltinsen30,  for 
example).  Under  potential  flow  theory  (assuming  inviscid,  irrotational  flow),  the  fluid 
velocity  within  the  fluid  domain  may  be  represented  as  the  gradient  of  a  scalar  potential 
function 

__,     .  cO      .  <?0>     .  c<D 

v  =  VcD  =  i^  +  j  —  +  k^^  (2-47) 

ex        dy         dz 

where  v  is  the  fluid  velocity  vector,  0(x,y,z,t)  is  the  scalar  velocity  potential  function 

(written  here  in  terms  of  Cartesian  coordinates),  and  V  is  the  vector  differential  operator 

The  velocity  potential  must  satisfy  Laplace's  equation 

c20      c2<&     d<$> 

V;O^  +  ^  +  ^-  =  0  (2-48) 

dx~       dy~       dz' 
within  the  fluid  domain,  and  satisfy  boundary  conditions  specific  to  the  geometry  of  the 
application.  For  the  case  of  a  body  of  general  shape  moving  near  a  free  surface  in  a  fluid 
of  finite  depth  (figure  2-15),  the  velocity  potential  must  satisfy  the  following  boundary 
conditions: 

=  V  •  n    (on  the  surface  of  the  body)  (2-49a) 

dn 


28Newman.  op.  cit. 

29Chakrabarti.  S.K..  Hydrodynamics  of  Offshore  Structures.  Springer- Verlag.  London.  1987. 

30Faltinsen.  O.M..  Sea  Loads  on  Ships  and  Offshore  Structures,  Cambridge  University  Press.  UK,  1990. 
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«30 

dy 


0    (at  y  =  -h  (for  finite  water  depth)) 


(2-4%) 


c:0        r<J>  rO      1  ,  , 

— +  g~  +  2VcD-V —  +  -V<D-V(VO-Vd>   =0 
cV  c>  ct      2 

(at  y  =  n  (on  the  free  surface))  (2-49c) 

where  n  is  the  normal  vector  of  the  body  surface  (pointing  from  the  body  into  the  fluid 

domain),  V  is  the  velocity  vector  of  the  body,  g  is  the  acceleration  due  to  gravity,  r|  is  the 

free  surface  elevation  above  the  mean,  and  h  is  the  mean  depth  of  the  water    The  free 

surface  boundary  condition,  equation  (2-49c),  represents  the  complete  nonlinear  free 

surface  boundary  condition31.   The  linearized  free  surface  boundary  condition  (linearized 

about  the  mean  free  surface)  would  be  written: 


-^  +  g— =  0    (aty  =  0) 
cV  cy 


(2-49d) 


-y 


7 7 7 7 7 7 7 7 7        7        7        7" 


Figure  2-15.  Body  of  general  shape  moving  in  fluid  with  free  surface. 


3 'Newman,  op.  cit..  p.  247. 
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Additionally,  a  radiation  condition  must  be  established  (as  a  boundary  condition  limit  at 
infinity)  in  order  to  ensure  that  radiated  waves  propagate  away  from  the  body    For 
solution  of  equation  (2-48)  in  the  time  domain,  initial  conditions  are  also  required: 

4>->0 

>0    (as  t -»-*>)  (2-50) 

dt 

Solution  of  Laplace's  equation  (2-48)  for  the  velocity  potential,  O,  is  obtained  by  the 

method  of  integral  equations  using  Green's  theorem  which  can  be  written  in  the  general 

form: 

J(cDV:G-GV:0)dV  =  J(<DVG-GVO)-ndS  (2-51) 

V  s 

where  V  is  defined  as  the  volume  bound  by  surface  S,  where  S  includes  the  body  surface, 
the  free  surface,  the  bottom  surface,  and  the  surface  at  "infinity"  (figure  2-16)    An 
appropriate  Green  function  (G)  can  be  obtained  which  satisfactorily  represents  the  source 
potentials  and  their  derivatives  for  various  boundary  condition  combinations  (Wehausen 
and  Laitone32,  Newman33,  for  example). 

Once  the  velocity  potential  function  is  obtained,  the  fluid  velocity  vector  in  the 
fluid  domain  can  be  found  by  definition  of  the  potential  function  as  written  in  equation 
(2-47).  The  hydrodynamic  pressure  exerted  on  the  body  by  the  fluid  is  obtained  from  the 
unsteady  Bernoulli  equations: 


TO      1  i 

— +  - VO 

at    21 


(2-52) 


where  P  is  the  hydrodynamic  pressure  (a  function  of  position)  and  p  is  the  fluid  density. 
The  hydrodynamic  force  exerted  on  the  body  (or  the  panel  of  interest)  by  the  fluid  is 
obtained  by  integrating  the  hydrodynamic  pressure  over  the  submerged  portion  of  the 
surface: 


F  =  JJ(P-n)dS  (2-53) 


where  F  is  the  force  vector  on  the  surface  or  panel,  S  is  the  submerged  surface  of  the  body 
or  panel,  and  n  is  the  direction  normal  of  the  panel  surface. 


32Wehausen.  J.V..  and  Laitone.  E.V..  "Surface  Waves".  Encyclopedia  of  Physics  (Vol.  9).  Spnnger- 
Verlag.  Berlin.  1960. 

33Ne\vman.  J.N..  "Algorithms  for  the  free  surface  Green  function".  Journal  of  Engineering  Mathematics. 
19.  1985.  pp.  57-67. 
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Figure  2- 1 6.   Surfaces  of  integration  for  Green's  theorem 


For  slender  bodies  (where  length  is  much  greater  than  width),  application  of 
two-dimensional  diffraction/radiation  potential  theory  has  proven  accurate  and  useful, 
particularly  at  higher  frequencies34  Several  two-dimensional  approaches  are  available  in 
the  literature,  the  most  common  of  which  is  strip  theory*-    In  strip  theory,  the  slender 
body  is  divided  into  many  thin  slices  or  strips,  each  of  which  has  approximately  constant 
cross-section,  so  that  two-dimensional  potential  theory  could  be  applied  to  determine  the 
approximate  hydrodynamic  forces  on  each  strip  independently  (figure  2-17).  For  rigid 
body  hydrodynamics  applications,  the  hydrodynamic  forces  on  each  strip  could  then  be 
integrated  along  the  length  of  the  body  to  determine  the  total  hydrodynamic  force  on  the 
body.  Alternatively,  similar  to  panel  methods,  the  hydrodynamic  forces  on  each  "strip" 
can  be  resolved  and  integrated  into  the  fluid-structure  interaction  dynamic  problem  as 
external  loads. 

Thus,  the  motions  of  a  large  structure  on  or  near  a  free  surface  (where  viscosity 
effects  can  be  ignored)  could  be  obtained  using  the  complete  potential  flow  theory 


34Lewis,  E.V.  (ed.).  Principles  of  \'a\'al  Architecture  (Vol.  Ill  -  Motions  in  li'ax'es  and  Controllability), 
Second  Revision.  SNAME.  Jersey  City.  NJ.  1989. 

35Salveson.  N..  Tuck.  E.O..  and  Faltinsen.  O..  "Ship  Motions  and  Sea  Loads".  SNAME  Transactions,  Vol. 
78.  1970.  pp.  250-279. 
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(linearized  or  fully  nonlinear)    It  computes  the  forces  exerted  on  the  structure  by  the 
incident  pressure  (the  external  "loading"  force  -  the  bubble  pulse  for  a  whipping  ship 
subjected  to  an  underwater  explosion  bubble),  the  diffraction  pressure  (scattering  of  the 
incident  pressure  from  the  "rigid"  structure),  and  the  radiation  pressure  (wave  radiation 
from  the  structure  as  it  moves  in  each  of  its  degrees  of  freedom)    Thus,  the  total  velocity 
potential  (<J>)  at  any  point  in  the  fluid  in  the  presence  of  a  moving  structure  (either  three- 
dimensional  if  panel  methods  are  used,  or  two-dimensional  if  strip  theory  is  used)  can  be 
written  as  a  superposition  of  the  velocity  potentials  due  to  incident,  diffracted,  and 
radiated  pressures  as  follows: 


\1       N 


d^o'+x^  +  ZE0:, 


(2-54) 


=1  1=1 


where  <I>0  is  the  velocity  potential  due  to  the  incident  pressure,  O*  is  the  velocity  potential 
due  to  the  diffracted  or  scattered  pressure  from  the  /th  panel  (where  M  is  the  total  number 
of  panels  making  up  the  total  surface  S),  and  3>[,  is  the  velocity  potential  due  to  the 

radiated  pressure  caused  by  motion  of  the  /th  panel  in  theyth  direction  (where  N  is  the 
number  of  degrees  of  freedom  for  each  panel)  . 


Figure  2-17.   Two-dimensional  strip  theory 
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For  most  marine  applications,  rigid  body  translations  are  ignored  and  the  solution 
of  the  equations  of  motion  is  carried  out  in  the  frequency  domain,  assuming  harmonic 
motion  of  the  structure  centered  around  an  equilibrium  position    The  incident  and 
diffracted  velocity  potentials  are  computed  on  the  assumption  that  the  structure  is  fixed  in 
its  equilibrium  position,  subjected  only  to  the  given  external  loading  (usually  wave 
loading)    The  radiated  velocity  potentials  are  then  computed  assuming  that  the  structure 
undergoes  one  simple  harmonic  motion  at  a  time  (i.e.  at  each  frequency  in  each  degree  of 
freedom)  in  calm  water  (see  Chakrabarti36,  Faltinsen37,  for  example)    For  assumed  simple 
harmonic  motion  of  the  form 

x  =  Xe"ot  (2-55) 

where  x  is  the  displacement  vector,  X  is  the  displacement  magnitude  vector,  and  co  is  the 
frequency,  the  total  hydrodynamic  force  on  the  structure  (calculated  with  potential  flow 
theory)  has  the  form: 

Fn  =  -»^„Mm  +  icoXC  (2-56) 

n  n         mn  n       mn  y  ' 

in  which  Fn  is  the  force  on  the  structure  in  the  wth  degree  of  freedom,  and  X    is  the 
displacement  of  the  structure  in  the  nth  degree  of  freedom  M^  is  called  the  added  mass 
matrix  (hydrodynamic  added  mass  in  the  nth  degree  of  freedom  due  to  unit  displacement 
in  the  mth  degree  of  freedom),  and  thus  the  first  term  on  the  right  hand  side  represents  the 
hydrodynamic  forces  which  are  in  phase  with  acceleration    C^  is  called  the  radiation 
damping  matrix  (hydrodynamic  damping  in  the  nth  degree  of  freedom  due  to  unit 
displacement  in  the  mth  degree  of  freedom),  and  thus  the  second  term  on  the  right  hand 
side  represents  the  hydrodynamic  forces  which  are  in  phase  with  velocity. 

For  applications  where  simple  harmonic  motion  about  a  fixed  equilibrium  position 
is  not  a  reasonable  assumption  (e.g.  when  rigid  body  translations  are  important),  the 
separate  accounting  for  added  mass  and  radiation  damping  forces  as  above  is  normally  not 
possible  (except  in  the  case  of  steady  rigid  body  translation).  For  cases  with  unsteady 
rigid  body  translation  in  addition  to  harmonic  vibration  (e.g.  whipping  of  a  submerged 
submarine  subjected  to  an  underwater  explosion  bubble),  the  calculation  of  the  total 
hydrodynamic  forces  must  be  accomplished  in  the  time-domain.  Application  of  potential 
flow  theory  in  the  time-domain  to  predict  general  motions  of  structures  is  not  new  and 
fundamentally  straight  forward,  but  has  had  limited  application  for  complex  geometries 
and  loadings  because  of  limited  computational  resources    Discussions  of  direct  time- 


36Chakrabarti.  op.  at. 
37Faltinsen.  op.  at. 


41 


domain  solution  techniques  are  presented  by  various  authors  for  various  applications 
(Beck  and  Magee38,  Ogilvie39,  for  example). 
2.4.2     Viscous  fluid  damping. 

As  mentioned  previously,  dynamic  energy  can  be  lost  into  the  surrounding  fluid 
medium  through  forces  associated  with  fluid  viscosity    First,  viscous  shear  forces  occur 
within  the  viscous  boundary  layer  surrounding  the  moving  structure  (so-called  "skin 
friction"  drag)    Second,  separation  of  the  boundary  layer  from  the  structure  leads  to 
additional  forces  on  the  moving  structure  (so-called  "form"  drag).  Together,  boundary 
layer  shear  forces  and  boundary  layer  separation  forces  constitute  the  viscous  fluid  forces 

The  extent  to  which  each  of  these  forces  is  important  is  a  function  of  the  geometry 
of  the  body,  and  the  relative  fluid  motion  around  the  body  (characteristic  amplitude  of 
relative  fluid  motion  and  characteristic  frequency)  Where  geometry  and  fluid  flow  around 
the  structure  are  such  that  the  boundary  layer  remains  attached  (eg  streamlined  bodies  at 
low  angles  of  attack),  boundary  layer  shear  forces  tend  to  dominate    However,  for  "bluff' 
bodies  such  as  circular  cylinders  and  "flat"  plates,  boundary  layer  separation  readily  occurs 
which  tends  to  produce  forces  greatly  exceeding  the  boundary  layer  shear  forces. 
Numerous  experiments  have  been  conducted  on  various  geometries  giving  indication 
under  which  conditions  each  of  the  forces  would  be  significant,  and  estimation  of  the 
forces. 

The  largest  body  of  data  exists  for  harmonic  transverse  oscillation  of  horizontal 
cylinders  for  which  two  dimensional  flow  around  the  body  can  be  assumed.  As  discussed 
previously,  flow  around  a  two-dimensional  body  is  typically  characterized  by  several  non- 
dimensional  parameters,  including  the  Keulegan-Carpenter  number  and  the  Reynolds' 
number  (or  alternatively  the  viscous-frequency  parameter,  3).  Anaturk40  and  Anaturk  et. 
al.41  present  regions  in  the  K.C  -  P  plane  where  two  dimensional  flows  possess  different 
characteristics    For  example,  when  KC  <  1  and  p  <  103  -  104  (Rn  <  103  -  104),  the 
boundary  layer  flow  typically  remains  attached  for  the  case  of  a  smooth  circular  cylinder 


38Beck.  R.F.,  and  Magee.  A.R..  "Time-domain  Analysis  for  Predicting  Ship  Motions".  Dynamics  of 

Marine  I'ehicles  and  Structures  in  Wa\'es.  Elsevier  Science  Publishers.  B.V..  1991.  pp  49-64. 

390gilvie.  T.F.,  "Recent  Progress  Toward  the  Understanding  and  Prediction  of  Ship  Motions". 

Proceedings  of  the  Fifth  Symposium  on  Na\'al  Hydrodynamics.  Office  of  Naval  Research.  Washington. 

DC.  1964.  pp.  3-128. 

40 Anaturk.  A.R..  "An  experimental  investigation  to  measure  hydrodynamic  forces  at  small  amplitudes  and 

high  frequencies".  Applied  Ocean  Research.  13.  4,  1991.  pp.  200-208. 

41  Anaturk.  A.R..  Tromans.  PS.,  van  Hazendonk.  H.C..  Sluis.  CM.,  and  Otter.  A..  "Drag  forces  on 

cylinders  oscillating  at  small  amplitude:  a  new  model".  Journal  of  Offshore  Mechanics  and  Arctic 

Engineering.  114.  1992.  pp.  91-103. 
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As  KC  and  (3  (or  Rn)  increase,  the  turbulence  in  the  boundary  layer  increases  and  flow 
separates 

For  the  case  where  the  two  dimensional  flow  does  not  separate,  the  flow  (and 
resulting  forces)  can  be  satisfactorily  predicted  using  an  attached  laminar  boundary  layer 
theory  (Newman42,  Stuart43,  Batchelor44,  for  example)    This  theory  predicts  the  drag 
forces  resulting  from  the  combined  normal  and  tangential  stresses  in  the  fluid  using 
Navier-Stokes  equations. 

For  a  body  of  general  shape,  the  forces  associated  with  separated  flow,  and  the 
point  at  which  transition  to  separated  flow  occurs  (i.e.  the  stability  of  the  boundary  layer) 
is  dependent  upon  many  factors    Many  attempts  have  been  made  to  solve  separated  flow 
around  marine  structures  numerically  (generally  limited  to  two  dimensional  flow). 
Examples  of  various  "state-of-the-art"  methods  include: 

(a)  Single  vortex  method  (Faltinsen  and  Sortland45) 

(b)  Vortex  sheet  model  (Faltinsen  and  Pettersen46) 

(c)  Discrete  vortex  method  (Sarpkaya47,  Clements48) 

(d)  Vortex-in-cell  method  (Stansby  and  Dixon49) 

(e)  Combination  of  Chorin's  method  and  vortex-in-cell  method  (Chorin50, 
Smith  and  Stansby51) 

(f)  Navier-Stokes  solvers  (Lecointe  and  Piquet52) 

Additionally,  attempts  have  been  made  to  include  vortex  shedding  models  in  numerical 
solutions  for  potential  flows  and  thus  combine  wave  radiation  with  the  effects  of  viscosity 


42Ne\vman.  op.  cit. 

43Stuart.  J.T..  "Double  boundary  layers  in  oscillatory  viscous  flow".  Journal  of  Fluid  Mechanics,  24. 

1966.  pp.  673-687. 

44Batchelor.  G.K.,  An  Introduction  to  Fluid  Dynamics.  Cambridge  University  Press.  London.  1973. 

45FalUnsen.  O.M  .  and  Sortland.  B..  "Slow-drift  eddy  making  damping  of  a  ship".  Applied  Ocean 

Research,  9.  1.  1987.  pp.  37-46. 

46Faltinsen.  O.M..  and  Pettersen.  B..  "Application  of  a  vortex  tracking  method  to  separated  flow  around 

marine  structures".  Journal  of  Fluids  and  Structures.  1.  1987.  pp.  217-237. 

47Sarpkaya.  T..  "Computational  Methods  With  Vortices  -  The  1988  Freeman  Scholar  Lecture".  Journal  of 

Fluids  Engineering.  Vol.  111.  March  1989.  pp.  5-45. 

48Clements.  R.R..  and  Maull.  D.J..  "The  representation  of  sheets  of  vorticity  by  discrete  vortices". 

Progress  in  Aerospace  Science.  Vol.  16.  No.  2.  1975,  pp.  129-146. 

49Stansby.  P.K..  and  Dixon.  AG..  "Simulation  of  flows  around  cylinders  by  a  Lagrangian  vortex  scheme". 

Applied  Ocean  Research.  Vol.  5.  No.  3.  1983.  pp.  167-178. 

50Chonn.  A.  J..  "Numerical  study  of  slightly  viscous  flow" '.  Journal  oj  Fluid Mechanics.  57.  1973. 

pp.  785-796. 

51Smith.  PA.,  and  Stansby.  P.K..  "Impulsively  started  flow  around  a  circular  cylinder  by  the  vortex 

method".  Journal  of  Fluid  Mechanics.  194.  1988.  pp.  45-77. 

52Lecointe.  Y..  and  Piquet.  J  "Compact  finite-difference  methods  for  solving  incompressible 

Navier-Stokes  equations  around  oscillatory  bodies".  I  'on  Karman  Institute  for  Fluid  Dynamics  Lecture 

Series  1985-04,  Computational  Fluid  Dynamics.  1985. 
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(Yeung  and  Wu5\  Wong  and  Calisal54,  for  example)    While  these  methods  have 
documented  satisfactory  agreement  with  experiment  in  some  cases,  they  are  generally  not 
considered  robust  enough  (numerical  instabilities,  etc  )  to  be  applied  with  complete 
confidence  to  new  problems  where  there  is  no  guidance  from  experiments55    Thus, 
numerical  methods  to  characterize  separated  flow  are  often  used  in  conjunction  with 
experimental  methods 

Morison's  equations    In  practice,  when  drag  forces  are  important  (and  wave 
radiation  may  be  ignored),  and  two-dimensional  oscillatory  flow  can  be  assumed,  a 
Morison's  equation  or  modified  Morison's  equation  can  be  used  to  predict  hydrodynamic 
forces  .  Originally  proposed  by  Morison  et  al.56  for  the  prediction  of  wave  forces  on 
vertical  piles  which  extend  from  the  bottom  through  the  free  surface,  the  Morison 
equation  has  been  modified  to  include  application  to  cylinders  of  varying  orientation  under 
various  loading  conditions,  including  cases  where  the  cylinders  are  free  to  move  under  the 
imposed  fluid  loading. 

Morison's  equation  for  a  stationary-  cylinder  in  an  oscillating  flow  field  includes 
the  effect  of  fluid  inertia  forces  and  assumes  drag  forces  to  be  a  quadratic  function  of 
velocity: 


PACmu  +  -pDCdu|u|  (2-57) 


where  F  is  the  hydrodynamic  force  (per  unit  length  of  cylinder  -  similar  to  "strip  theory", 
the  total  force  on  the  cylinder  can  be  found  by  integrating  over  the  length  of  the  cylinder), 
A  is  cross  sectional  area,  D  is  diameter,  u  is  fluid  velocity  (li  is  fluid  acceleration), 
coefficient  Cm  is  called  the  inertia  coefficient  (the  first  term  being  proportional  to 
acceleration  is  equivalent  to  an  inertia  force),  and  coefficient  Cd  is  called  the  drag 
coefficient  (the  second  term  being  proportional  to  velocity  (squared)  is  equivalent  to  a 
drag  force).  Figure  2-18  illustrates  the  parameters  used  in  Morison's  equation 

Morison's  equation  for  an  oscillating  cylinder  in  a  stationary- fluid  includes  the 
effect  of  fluid  added  mass  forces  and  again  assumes  drag  forces  to  be  a  quadratic  function 
of  velocity,  giving  the  reaction  forces  on  the  cylinder: 


53Yeung.  R.W..  and  Wu.  C.  "Viscosity  effects  on  the  radiation  hydrodynamics  of  two-dimensional 

cylinders".  Proceedings  of  the  10th  International  Conference  on  Offshore  Mechanics  and  Arctic 

Engineering.  ASME.  1991.  pp.  309-316. 

54Wong.  L.H..  and  Calisal.  S.M..  "A  numerical  solution  for  potential  flows  including  the  effects  of  vortex 

shedding".  Proceedings  of  the  1 1th  International  Conference  on  Offshore  Mechanics  and  Arctic 

Engineering.  ASME.  1992.  pp.  363-368. 

55Faltinsen.  op.  at. 

56Monson.  JR..  O'Brien.  MP.  Johnson.  J.W..  and  Schaaf.  S.A..  "The  force  exerted  by  surface  waves  on 

piles".  Petroleum  Transactions.  AIME.  189.  1950.  pp.  149-157. 
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-pACAx--pDCdx|x| 


(2-58) 


where  x  and  x  are  the  acceleration  and  velocity  of  the  cylinder,  CA  is  called  the  added 
mass  coefficient,  and  Cd  is  the  drag  coefficient  for  and  oscillating  cylinder  in  still  water 


Length,    L 


cylinder 
velocity,   x 
acceleration,  x 


"far  field" 
water  particle 
velocity,   acceleration 


P  fluid   density 

v   fluid   kinematic   viscosity 


Figure  2-18.  Parameters  used  in  Morison's  equation 

Coefficients  Cm,  Cd,  CA,  and  Cd  are  normally  experimentally  determined  for 
various  geometries  and  loading  conditions  (usually  as  functions  of  KC  and  (3  or  Rn),  and 
have  been  determined  in  a  large  number  of  experimental  studies  (Bearman  et.  al.57, 
Sarpkaya58,  Chakrabarti  and  Cotter59,  for  example). 

When  a  cylinder  is  free  to  move  in  a  moving  fluid,  Morison's  equation  can  be 
written  in  an  extended  form  by  combining  equations  (2-57)  and  (2-58): 


F  =  pACmu  +  —  pDCdu|u|-pACAx  — pDCdx|x| 


(2-59) 


2"        -   '   '     ■        "       2 
This  form  is  known  as  the  independent  flow  fields  model60  and  is  obtained  as  the 
superposition  of  the  two  independent  flow  fields,  a  far  field  due  to  the  water  particle 
motion  (i.e.  the  external  loading)  and  relatively  unaffected  by  the  presence  and  motion  of 


57Bearman.  P.W..  Lin.  X.W..  and  Mackwood.  PR..  "Measurement  and  prediction  of  response  of  circular 

cylinders  in  oscillating  flow".  Proceedings  of  the  6th  International  Conference  on  the  Beha\'ior  of 

Offshore  Structures  (BOSS).  1992.  pp.  297-307. 

58Sarpkaya.  T..  "Force  on  a  cylinder  in  viscous  oscillatory  flow  at  low  Keulegan-Carpenter  numbers". 

Journal  of  Fluid  Mechanics,  Vol.  165,  1986.  pp.  61-71. 

59Chakrabarti.  S.K..  and  Cotter.  DC.  "Hydrodynamic  coefficients  of  mooring  tower".  Proceedings  of  the 

15th  Annual  Offshore  Technology  Conference,  OTC  4496.  1983.  pp.  449-458. 

60Chakrabarti.  op.  at.,  p.  187. 
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the  structure,  and  a  near  field  resulting  from  the  motion  of  the  structure    Coefficients  Cm 
and  Cd  may  be  obtained  from  experiments  where  the  cylinder  is  held  fixed  in  oscillating 
fluid,  and  coefficients  CA,  and  Cj,  may  be  obtained  from  experiments  of  an  oscillating 
cylinder  in  otherwise  still  fluid    Alternatively,  the  coefficients  could  be  determined  from 
empirical  relations  (using  KC  and  Rn,  for  example): 

Cm,Cd=ffKCf=^,Rnf=-^l  (2-60) 


cA,ci 


(  x  T1  x  D^l 

1  KC    =i%L,Rn    =■ 


(2-61) 


D  n         v 

where  u0,x0  are  the  amplitudes  of  water  particle  and  structure  velocity,  respectively, 
T,T'  are  periods  of  oscillation  of  water  particle  and  structure,  respectively,  and  subscripts 
/and  n  refer  to  the  far  field  and  near  field,  respectively. 

An  alternative  to  the  independent  flow  fields  model  is  the  relative  velocity  model, 
in  which  forces  are  written  in  terms  of  the  relative  motion  between  the  structure  and  the 
fluid,  in  which  single  coefficients  are  assumed  to  apply: 

F  =  pACm  (u  -  x)  +  pAx  +  -pDCd  (u  -  x)|u  -  x|  (2-62) 

In  the  relative  velocity  model,  the  coefficients  may  be  determined  using  empirical  relations 
(using  KC  and  Rn,  for  example): 


(  v    T  v    D\ 

Cm,Cd=f   KCr=-^,Rnr- 


(2-63) 


D  '  v 

where  subscript  r  refers  to  the  relative  motion,  vr  is  the  relative  velocity  between  the 
structure  and  the  fluid  (vr0  is  the  amplitude  of  the  relative  velocity),  and  Tr  is  the 
combined  period  of  the  relative  motion    The  relative  velocity  model  has  seen  extensive 
use  in  the  analysis  of  hydrodynamic  forces  on  drag  dominated  offshore  structures  (Laya  et 
al.61,  Spidsoe  and  Karunakaran62,  for  example). 

Thus,  the  Morison  equation  can  be  used  under  various  conditions  to  predict  the 
hydrodynamic  forces  on  structures  (or  portions  of  structures)  where  drag  forces  and  fluid 


61Laya.  E.J..  Conner,  M..  and  Sunder.  S.S.,  "Hydrodynamic  forces  on  flexible  offshore  structures". 
Journal  of  Engineering  Mechanics,  ASCE.  March  1984.  pp.  433-448. 

62Spidsoe.  N..  and  Karunkaran.  D  .  "Nonlinear  effects  of  damping  to  dynamic  amplification  factors  for 
drag-dominated  offshore  platforms".  Proceedings  of  the  I  lth  International  Conference  on  Offshore 
Mechanics  and  Arctic  Engineering,  ASME.  1992.  pp.  231-239. 
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inertia  forces  are  .mportant    An  application  of  Monson's  equat.on  to  the  development  of  a 
model  for  submarine  whipping  will  be  made  in  chapter  3 
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3.  MODEL  FOR  SUBMARINE  WHIPPING  INCLUDING  DAMPING  EFFECTS 

3.1        Introduction. 

In  order  to  gain  some  quantitative  insight  into  the  relative  amount  of  damping 
which  could  be  provided  by  each  of  the  mechanisms  discussed  in  the  previous  chapter,  it  is 
first  necessary  to  develop  a  general  model  for  submarine  whipping,  which  may  include  any 
or  all  of  the  specific  mechanisms    One  way  of  doing  this  would  be  to  develop  a  simplified 
model  of  each  of  the  mechanisms  based  upon  the  fundamental  principles  of  mechanics, 
some  of  which  were  used  in  the  discussion  in  the  previous  chapter,  and  utilize  a  flexible 
formulation  for  the  equations  of  motion  which  permits  inclusion  of  selected  mechanisms 
Using  this  approach,  any  individual  mechanism  could  be  evaluated  separately  for  its 
expected  effect  on  total  damping 

The  equations  of  motion.  The  dynamic  behavior  of  any  structure  is  described  by 
the  equations  of  motion,  which  are  simply  extensions  of  Newton's  second  law  of  motion 
Fundamentally,  an  equation  of  motion  can  be  written  as  a  sum  of  forces  which  resist 
motion  of  the  structure  (i.e.  the  inertial  forces  (Finertldl ),  the  damping  or  dissipative  forces 

(^damping)'  anc*  tne  elastic  or  stiffness  forces  (Felastic )),  and  forces  which  induce  motion  (i.e. 

the  external  or  loading  forces  (Floading )).  Thus,  an  equation  of  motion  may  be  written  in  the 

general  form: 

F         +  F  +  F         =  F  (3-1  \ 

inertial  damping  elastic  loading  '  / 

Inertial  forces  are  those  forces  through  which  kinetic  energy  is  stored,  and  are 
dependent  upon  acceleration  and  a  mass  coefficient  ( Finertial  (M,  x) ).  Elastic  forces  are 
those  forces  through  which  potential  energy  is  stored  and  are  generally  dependent  upon 
displacement  and  some  elastic  or  stiffness  coefficient  (Felastie(K,x)).  As  discussed 
previously,  damping  forces  are  those  forces  related  to  the  dissipation  of  energy,  and 
depend  upon  the  mechanisms  through  which  energy  is  dissipated,  generally  dependent 
upon  velocity  and  some  damping  coefficient  (Fdampmg  (C,  x)).    It  should  be  noted  that,  in 

general,  these  dependencies  are  neither  constant  nor  linear 

The  external  loading  forces  include  all  forces  which  are  "externally"  applied  to  the 
structure.  These  forces  may  themselves  be  thought  of  as  being  inertial,  dissipative,  or 
elastic  in  nature    Thus,  equation  (3-1 )  may  be  rewritten  by  redefining  the  inertial, 
damping,  and  elastic  terms  such  that  they  also  include  the  "external"  loading  forces.  For 
example,  for  a  system  consisting  of  a  structure  moving  in  a  fluid,  equation  (3-1 )  may  be 
rewritten 
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F1        +  F'  +  F1        -  ft  n.  ~>\ 

1  inertial  *  damping      '    L  elastic    —  u  {->-£) 

where  F,L,ai  represents  all  inertial  forces  (inertia  of  the  structure  and  inertia  of  the 
surrounding  fluid),  Fj^^  represents  all  damping  or  dissipative  forces  (energy  dissipation 
within  the  structure  and  energy  lost  into  the  surrounding  fluid),  and  Fe'lastic  represents  all 
elastic  forces  (structural  stiffness  and  hydrodynamic  buoyancy  forces) 

In  general,  the  terms  associated  with  inertial  forces,  damping  forces,  and  elastic 
forces  may  be  written  in  any  order  as  long  as  all  forces  are  accounted  for  (and  as  long  as  a 
proper  sign  convention  is  maintained)    In  keeping  with  the  classification  of  damping 
forces  discussed  in  chapter  2,  the  equations  of  motion  for  a  whipping  ship  might  therefore 
be  written  in  the  form: 

''hull  +  ^internal   +  ^external    =  ^  (->-->) 

where: 

f=(f       +  f        +  f      ) 

hull  I     inertial  damping  elastic  J  ... 

C  -fp  +p  +p  ] 

int  emal  I     inertial  damping  elastic  J  . 

p  -  \  p  +  p  +  p  I  CX-A\ 

external  I     inertial  damping  elastic  J  .  \J   ~ ) 

where  Fe>;terTial  would  include  all  external  forces  on  the  hull.  Thus,  the  dynamics  of  a 
whipping  ship  may  be  thought  of  as  a  (coupled)  superposition  of  the  dynamics  of  the  hull, 
any  internal  components,  and  the  external  fluid  medium    Equations  (3-3)  and  (3-4) 
represent  one  form  of  the  complete,  coupled,  nonlinear  equations  of  motion  for  a 
whipping  ship. 

Using  this  approach  to  the  equations  of  motion,  each  term  may  be  represented  by 
an  independent  formulation  (which  may  involve  specific  assumptions  of  linearity  for  that 
term),  so  long  as  the  terms  remain  coupled  via  a  single  solution  technique    As  an  example, 
the  dynamic  behavior  of  a  slender  ship  hull  (such  as  a  modern  submarine)  might  be 
calculated  based  upon  a  linear  simple  beam  finite  element  model,  while  the  external 
(hydrodynamic)  forces  might  be  calculated  based  upon  a  nonlinear  three-dimensional  panel 
method.  As  will  be  discussed  subsequently,  even  if  the  hull  itself  can  be  assumed  to 
respond  linearly,  the  nonlinearities  associated  with  the  external  (hydrodynamic)  forces 
require  the  solution  of  the  equations  of  motion  in  the  time-domain.  It  is  the  existence  of 
the  coupling  of  the  terms  within  the  equations  of  motion  which  require  the  solution  of  the 
combined  equation  by  a  single  numerical  algorithm. 
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3.2        Formulation  of  the  dynamic  analysis  procedure. 

In  general  structural  dynamics,  the  selection  of  the  most  appropriate  method  to 
solve  the  equations  of  motion  depends  upon  the  type  of  loading  on  the  structure,  the 
extent  to  which  the  equations  of  motion  can  be  uncoupled  and  linearized,  and  the 
existence  of  rigid  body  motions    Nonlinearities  which  must  be  addressed  include 
geometric  nonlinearities  (large  rotations),  material  nonlinearities  (nonlinear  material  stress- 
strain  relations),  nonlinearities  associated  with  the  external  loading  (hydrodynamic 
nonlinearities),  nonlinear  boundary'  conditions,  and  so  on.  The  extent  to  which 
nonlinearities  can  be  "linearized"  depends  upon  the  assumptions  made  in  developing  the 
physical  model  for  the  motions  of  the  particular  structure    For  ship  whipping,  even  if 
material  and  geometric  nonlinearities  can  be  effectively  linearized  (by  assuming  linear 
elastic  stress-strain  relations  and  small  rotations  -  reasonable  assumptions  for  the  case  of 
elastic  whipping  of  a  steel  ship  hull),  hydrodynamic  forces  can  be  expected  to  be  highly 
non-linear  (as  discussed  in  chapter  2),  and  thus  complete  linearization  of  the  equations  of 
motion  is  generally  not  possible    Additionally,  for  the  general  dynamic  response  of  a 
whipping  ship  subjected  to  an  underwater  explosion  bubble,  it  could  be  expected  that  there 
would  be  a  significant  potential  for  rigid  body  motion    For  these  reasons,  the  equations  of 
motion  must  be  solved  directly  in  the  time-domain  (utilizing  a  time-stepping  numerical 
algorithm). 

Implementation  of  time-domain  analysis  procedures  for  the  dynamic  response  of  a 
structure  in  a  fluid  medium  is  discussed  in  great  detail  in  a  large  number  of  texts  and 
journals  articles.  Bathe63,  Clough  and  Penzien64,  and  Argyris  and  Mlejnek65  provide 
particularly  thorough  discussions  of  the  solution  of  nonlinear  equations  of  motion  by  direct 
time  integration  methods    Direct  time  integration  means  that  the  equations  of  motion 
would  be  solved  using  a  numerical  step-by-step  procedure.  The  term  "direct"  generally 
implies  that  there  is  no  transformation  of  the  equations  into  a  different  form  prior  to 
numerical  integration,  although  using  the  formulation  of  the  equations  of  motion  discussed 
above,  it  is  feasible  to  utilize  a  coordinate  transformation,  so  long  as  it  is  conducted  at 
each  time  step  (this  will  be  discussed  in  greater  detail  later).  In  essence,  instead  of  trying 
to  satisfy  equilibrium  at  only  one  distinct  time  /  (as  for  static  equilibrium),  direct  time 
integration  aims  to  satisfy  dynamic  equilibrium  at  all  discrete  time  intervals  At  apart    This 
implies  that,  basically,  dynamic  equilibrium  (which  includes  the  effects  of  inertia,  damping 


63Bathe.  op.  at. 

64Clough  and  Penzien.  op.  at. 

65Argvns.  J.,  and  Mlejnek.  J. P..  Dynamics  of  Structures.  Elsevier  Science  Publishers.  The  Netherlands. 
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and  elastic  forces)  is  sought  at  discrete  time  points  within  the  interval  of  solution 
Typically,  for  direct  time  integration  algorithms,  the  solution  at  the  "next"  required  time, 
/+A/,  is  calculated  based  upon  the  solutions  at  all  previous  times,  and  assuming  an 
incremental  variation  in  displacements,  velocities,  and  accelerations    It  is  the  form  of  the 
assumption  on  the  variation  of  displacements,  velocities,  and  accelerations  within  each 
time  interval  that  determines  the  accuracy,  stability,  and  required  time  of  the  solution 
procedure. 

For  direct  time  integration  methods,  the  incremental  equation  of  motion  can  be 
written  as  an  extension  of  equation  (3-2)  by  subtracting  the  equilibrium  equation  at  time  t 
from  that  at  time  t+At: 

[p.fAt       _pt  ),fpt+it  _pt  1  (pt+At     _pl  ]   =   0  r,« 

I     inertial  inertial  J  I     damping  damping  J  I     elastic  elastic  J  \~    J ) 

As  discussed  above,  inertial  forces  are  dependent  upon  acceleration  and  a  mass  coefficient 
(Finertial  (M,  x)),  damping  forces  are  dependent  upon  velocity  and  some  damping  coefficient 

(Fdamping(C,x)),  and  elastic  forces  are  dependent  upon  displacement  and  some  elastic  or 
stiffness  coefficient  (Felastlc(K,x))    Thus,  if  coefficients  M,  C,  and  K  can  be  related  to 
acceleration,  velocity  and  displacement  (respectively),  then  equation  (3-5)  contains  only 
three  unknowns  at  time  t+At  (the  solution  at  time  t  would  be  known)  and  may  be  solved  if 
two  of  the  three  are  known  as  relations  to  the  third 

The  Newmark  time  integration  method.  The  Newmark  integration  method66  , 
which  is  one  of  the  most  reliable  direct  time  integration  schemes,  and  whose  stability 
condition  is  well  known,  can  be  employed  quite  appropriately  for  this  application  (see  also 
Bathe67,  Argyris  and  Mlejnek68,  and  Doyle69)    In  the  Newmark  method  (often  referred  to 
as  the  Newmark  (3  Method),  iterative  calculation  is  carried  out  at  each  time  step  and  is 
moved  to  the  next  time  step  after  an  equilibrium  condition  is  reached. 

The  basic  equations  of  the  Newmark  method  are  derived  by  assuming  Taylor 
expansions  for  displacements  (x)  and  velocities  (x ): 

xt+,t=xt+xtAt  +  xt^-  +  P(xt+At-xt)At2  (3-6) 

xt+At=xt+xtAt  +  y(xt+At-xt)At  (3-7) 


66Ne\vmark.  N.M..  "A  method  of  computation  for  structural  dynamics."  American  Society  of  Civil 

Engineers  Transactions.  Vol.  127.  1962.  pp.  601-630. 

67Bathe.  op.  at. 

68 Argyris  and  Mlejnek.  op.  tit.,  pp.  485-488. 

69Doyle.  JR.  Static  and  Dynamic  Analysis  of  Structures.  Kluwer  Academic  Publishers.  Boston.  1991.  pp. 

353-357. 
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where  (3  and  y  are  coefficients  to  evaluate  contributions  of  remainders.  Thus,  the 
displacements  and  velocities  are  approximated  based  upon  "assumed"  accelerations  at 
each  time  interval  (thus  requiring  an  iterative  approach),  and  based  upon  the  solution  at 
the  previous  time  interval.  Newmark  showed  that  y  and  (3  are  not  entirely  arbitrary  and 
their  selection  determines  the  stability  and  accuracy  of  the  numerical  solution    Newmark 
found  that  unless  y  is  greater  than  or  equal  to  0.5,  spurious  damping  is  introduced  into  the 
equations  of  motion    Additionally,  unless  (3  is  greater  than  or  equal  to  0.25,  there  is  an 
incorrect  maximum  velocity  response    Further,  with  y  =  0  5  and  (3  =  0.25,  the  scheme 
leads  to  unconditionally  stable  solutions 

To  implement  the  Newmark  method  in  the  analysis  of  the  dynamics  of  a  structure 
subject  to  nonlinear  loading  forces  (such  as  a  whipping  submarine),  the  equations  of 
motion  are  first  put  in  the  form: 

resisting  loading  w"°7 

where  Fresi       refers  to  all  those  forces  associated  with  the  structure  itself  which  resist 

motion  (i.e.  Fhull  as  defined  in  section  3.1)  and  Floading  refers  to  all  (other)  forces  which 

induce  motion  (i.e.  F,nlemal  +  Fextemal  as  defined  in  section  3  1)    Defining  Fresistmg  for  a 

structure  as: 

F^ing  =  Mx  +  Cx -f  Kx  (3-9) 

the  equations  of  motion  for  a  structure  can  be  written  in  the  general  form: 

Mx  +  Cx  +  K.x  =  F,oadmg  (3-10) 

or  in  incremental  form: 

M(x\.  v  -  Xt  )  +  C(Xt+A,  ~Xt)+  K(Xt-  V  -  Xt  )  =  Floadmgt+At  -  Fload,ngt  C3'1  ]  ) 

From  equation  (3-10),  the  acceleration  required  for  dynamic  equilibrium  to  exist  at 
time  t+At  can  be  written: 

K*  =M-,{{FllwdlJ^t  -{Cxt+At  +  Kxt+At}}  (3-12) 

Because  Fload     is  generally  also  dependent  upon  the  accelerations  (and  displacements  and 

velocities  by  equations  (3-6)  and  (3-7)),  it  is  evident  that  an  iterative  procedure  must  be 
implemented  at  each  time  t+At  in  order  to  obtain  equilibrium  by  equation  (3-12). 
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Newmark70  proposed  a  general  method  for  the  solution  of  structural  dynamics 
which  may  be  summarized  as  follows: 

1.  Assume  initial  values  of  displacement,  velocity,  and  acceleration  (i  e  at 
each  node)  for  time  t  =  0. 

2.  For  each  time  step: 

(a)  Assume  values  of  acceleration  for  the  current  time  step 

(b)  Compute  displacement  and  velocity  for  the  current  acceleration 
(from  equations  (3-6)  and  (3-7). 

(c)  For  the  current  displacement,  velocity,  and  acceleration,  compute 
the  loading  and  resisting  forces  in  accordance  with  the  physical 
model  for  the  particular  application. 

(d)  Compute  the  acceleration  required  to  obtain  dynamic  equilibrium 
with  these  loading  and  resisting  forces  (from  equation  (3-1 1 )) 

(e)  Compare  the  required  acceleration  with  the  assumed  acceleration 
for  the  current  time  step.  If  these  are  the  same  (i.e.  within 
tolerance),  the  calculation  for  this  time  step  is  complete    If  these 
are  different,  repeat  the  calculation  (i.e.  perform  another  iteration 
through  steps  (a)  -  (e))  with  a  different  value  of  assumed 
acceleration  (the  derived  acceleration  could  be  used  as  the  new 
assumed  acceleration). 

3.  Move  onto  the  next  time  step  (t+At),  assuming  the  final  acceleration  for  the 
last  time  step  (t). 

It  should  be  briefly  noted  here  that  the  modeling  of  internal  damping  devices  such 
as  structures  or  liquids  requires  only  the  addition  of  TV  equations  of  motion  to  the  set 
defined  above  (for  an  internal  damping  device  with  TV  participating  modes  of  response). 
Thus,  in  addition  to  the  requirement  for  equations  (3-12)  to  be  satisfied  at  each  time  t+At, 
the  following  equations  must  also  be  satisfied  at  each  time  t+At: 

x    =-£n£iniL  (3-13) 

where  subscript  n  is  for  each  mode  of  the  internal  damping  device  (up  to  mode  N).  Thus, 
if  there  are  M  degrees  of  freedom  required  to  represent  the  displacements  of  the  hull,  and 
N  (modal)  degrees  of  freedom  required  to  represent  the  (modal)  displacements  of  the 
internal  device,  then  a  total  ofM-N  simultaneous  equations  must  be  satisfied  using  the 
form  of  equations  (3-12)  and  (3-13). 


70Ne\vmark.  op.  cit..  pp.  606-611. 
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To  implement  the  Newmark  method  in  the  analysis  of  the  dynamics  of  a  structure 
subject  to  loading  that  is  highly  nonlinear  (such  as  a  whipping  submarine  subject  to 
nonlinear  hydrodynamic  forces),  a  slight  modification  to  the  above  method  must  be  made 
Because  of  the  highly  nonlinear  loading  forces,  the  equations  of  motion  must  be  solved  at 
each  time  step  using  an  iterative  process  (such  as  a  Newton-Raphson  process)  which 
"guarantees"  convergence.  Without  such  an  iterative  process,  the  solution  may  not 
converge  and  equilibrium  may  not  be  (numerically)  obtained 

Suzuki  and  Yoshida71  discuss  application  of  the  Newmark  method  with  solution  at 
each  time  step  by  an  iterative  Newton-Raphson  method.   In  the  case  of  elastic  ship 
whipping  with  nonlinear  hydrodynamic  loading,  the  specific  formulation  of  the  Newton 
iteration  may  proceed  by  using  equations  (3-6),  (3-7)  and  (3-1 1),  and  expressing 
acceleration  and  velocity  at  time  t+At  in  terms  of  displacement  at  time  t+At 


1 


L  t+At 


PAt: 


(xt+At  -xt) 


1 


( 


PAt 


x,  + 


i--L 

2(3 


\ 


(3-14) 


Y 


't-  m 


pAt 


(tw 


2Pj 


x,At 


(3-15) 


In  order  for  equilibrium  to  exist  at  time  t+At,  equations  (3-1 1),  (3-14)  and  (3-15)  must  be 
simultaneously  satisfied    Residuals  (which  must  vanish  for  equilibrium  to  exist)  may  be 
defined  for  equations  (3-11),  (3-14)  and  (3-15): 
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The  Newton  iteration  is  carried  out  such  that  the  unsatisfied  equilibrium  at  iteration^  is  to 
be  satisfied  at  iterationyW    In  other  words,  all  residuals  at  iteration^  I  must 
simultaneously  vanish.  To  obtain  expressions  for  residuals  at  iteration  stepy-  /,  the 


71  Suzuki.  R.  and  Yoshida.  K..  "Three  dimensional  nonlinear  dynamic  analysis  method  of  underwater  line 
structure  and  its  validation".  Proceedings  of  the  I  Oth  International  Conference  on  Offshore  Mechanics 
and  Arctic  Engineering  (OMAE),  ASME.  1991.  pp.  87-94. 
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residuals  at  iteration  stepy  may  be  expanded  in  linearized  Taylor  expansion  (in  matrix 
form): 
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(3-19) 
where  8x  =  (xt+Atj+]  -xt+Atj),  8x  =  (xt+Atj+]  -  xt+At . )  and  5x  =  (xt+Atj+]  -xt+at.)  are 

displacement,  velocity,  and  acceleration  corrections  to  be  made  for  they    7th  iteration  at 
time  step  t+At.   An  expression  for  the  displacement  correction  (5x)  may  be  obtained  by 
substituting  the  last  two  equations  in  (3-19)  into  the  first,  giving: 
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(3-20) 
Solving  equation  (3-20)  for  5x  and  substituting  this  term  back  into  equation  (3-19)  gives 
expressions  for  improved  displacement,  velocity  and  acceleration  (i.e.  displacement, 
velocity  and  acceleration  for  iteration  stepy    /  based  upon  iteration  stepy) 
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The  Jacobian  matrices, 
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,  may  be  approximated  at  each  iteration  step  by 
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With  these  modifications,  the  calculation  of  dynamic  equilibrium  using  the 
Newmark  method  for  the  case  of  highly  nonlinear  loading  could  proceed  as  follows: 
1  Assume  initial  values  of  displacement,  velocity  and  acceleration  (i.e. 

vectors  of  nodal  values) 

2.  For  each  time  step 

(a)  Assume  initial  values  of  displacement,  velocity,  and  acceleration 
equal  to  those  of  the  last  time  step 

(b)  Calculate  the  loading  forces  for  the  current  displacement,  velocity, 
and  acceleration  using  the  chosen  physical  model 

(c)  Calculate  the  Jacobian  matrices  (numerically)  for  the  current 
loading  forces 

(d)  Calculate  the  residuals  using  equations  (3-16)  through  (3-1 8). 

(e)  Solve  for  the  displacement  correction  (5x)  using  equation  (3-20) 

(f)  Calculate  improved  values  of  displacement,  velocity,  and 
acceleration  using  equations  (3-21). 

(g)  Reiterate  through  steps  (b)  through  (g)  until  the  residuals  are  within 
a  specified  tolerance  (i.e.  a  vector  norm  tolerance). 

3.  Move  onto  the  next  time  step,  assuming  the  final  displacements,  velocities, 
and  accelerations  of  the  previous  time  step 
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3.3        Definition  of  dynamic  forces. 

3.3.1      Hull  structure  forces. 

As  discussed  in  section  3.1,  the  forces  associated  with  the  dynamics  of  the  hull 
structure  may  be  represented  by  an  independent  formulation  from  that  used  to  represent 
the  internal  and  external  forces    This  formulation  may  involve  specific  assumptions  of 
linearity  which  may  be  applied  only  to  the  dynamics  of  the  hull  structure 

3.3.1.1  Finite  element  beam  model. 

In  representing  the  whipping  response  of  most  conventional  ships  to  low  frequency 
bubble  pulse  loading,  it  has  been  found  that  a  simple  beam  model  of  the  hull  structure  has 
provided  good  results72    This  idealization  is  supported  by  recognizing  that  most 
conventional  ships  (particularly  submarines)  have  length/beam  ratios  greater  than  eight 
Thus,  the  low  frequency  vibrational  modes  occur  as  beam-like  flexural  motions,  which  are 
represented  well  by  simple  beam  theory.   It  should  be  noted,  however,  that  nonlinear  hull 
response  (i.e.  where  the  assumptions  of  linear  elastic  stress-strain  relations  and  small 
rotations  cannot  be  made),  makes  the  utility  of  the  simple  beam  theory  limited    For 
nonlinear  hull  response,  a  fully  nonlinear,  three-dimensional  finite  element  technique  would 
be  more  appropriate    Therefore,  for  this  analysis  procedure,  the  assumptions  of  linear 
elastic  stress-strain  relations  and  small  rotations  will  be  made  in  order  to  justify  the  use  of 
the  simple  beam  model 

Because  sectional  properties  vary  gradually  along  the  length  of  the  hull,  a  finite 
element  beam  model  could  be  used  to  represent  the  structural  properties  of  the  submarine 
(figure  3-1).    The  model  consists  of  a  specified  number  of  beam  elements  (equal  to  the 
number  of  hull  stations  specified  on  the  ship's  drawings).  Most  ships  have  19  stations 
specified,  and  thus  would  be  modeled  with  19  beam  elements.  A  three-dimensional  beam 
model  will  be  used  because  of  the  directionality  of  the  loading  of  the  underwater 
explosion,  as  well  as  the  lack  of  complete  axial  symmetry  of  submarines  (i.e.  properties 
differ  in  each  plane  of  symmetry). 

Most  methods  for  structural  finite  element  analysis  develop  matrices  to  represent 
the  stiffness  and  mass  properties  of  each  finite  element,  and  then  combine  element  matrices 
into  global  or  structure  matrices    Modeling  of  beams  using  a  finite  element  technique,  and 
assignment  of  the  properties  of  each  finite  element  is  discussed  in  numerous  sources 


72Hicks.  op.  cit.,  p.  395. 
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I 

2. 
3. 


(Smith73,  Bathe74,  Doyle75,  Hughes76,  for  example)    An  appropriate  finite  element  analysis 
procedure  which  could  be  used  to  approximate  the  hull  forces  associated  with  submarine 
hull  whipping  could  be  outlined  as  follows. 

Define  the  model  (coordinate  system,  nodes,  element  types,  etc  ) 
Determine  each  element's  stiffness  and  mass  matrices  (ke  ,me)  and 
transform  them  to  hull  (or  global)  coordinates 
Assemble  the  hull  or  global  stiffness  and  mass  matrices  ( K,  M )  by 
superposition  of  the  element  matrices 

Approximate  a  reasonable  hull  damping  matrix  (C)  by  applying  appropriate 
Rayleigh  coefficients,  as  discussed  in  sections  2.2.3  and  3  3  12 
(C  =  aM  +  pK). 

Combine  the  hull  mass,  damping,  and  stiffness  matrices  to  obtain  the  total 
hull  force  as  described  in  section  3  1  (Fhull  =  Mx  +  Cx  +  Kx) 
This  procedure  will  be  incorporated  into  an  overall  computational  algorithm  and  used  for 
analyses  in  the  next  chapter 
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Figure  3-1 .   Structural  model  for  submarine  hull  whipping. 

Figure  3-2  illustrates  a  two-dimensional,  two-node  beam  element  with  vertical 
translation  and  rotation  as  the  nodal  displacements.  It  should  be  noted  that  because  a 
submarine  should  be  classified  as  a  relatively  thick  beam  (i.e.  the  cross-sectional 
dimensions  are  not  small  compared  to  the  length),  the  effects  of  rotary  inertia  and  shear 
deformation  at  each  element  should  be  accounted  for  in  the  finite  element  model 


73 Smith,  op.  cit. 
74Bathe.  op.  cit. 
75Doyle.  op.  cit. 
76Hughes.  OF..  Ship  Structural  Design.  SNAME.  Jersey  City.  NJ.  1988. 
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To  be  rigorous,  a  three-dimensional  beam  element  could  be  used  which  also  considers 
additional  nodal  displacements  in  the  horizontal  plane  as  well  as  torsion  and  axial 
displacements  (a  total  of  6  degrees  of  freedom  at  each  node).  Figure  3-3  illustrates  a 
general  three-dimensional,  two-node  beam  finite  element  with  6  degrees  of  freedom  at 
each  node.  It  should  be  noted  that  axial  displacements  can  generally  be  ignored  for 
submarine  whipping  due  to  the  lack  of  axial  loading  provided  by  the  surrounding  fluid. 
Additionally,  torsional  displacements  may  also  be  ignored,  although  it  could  be  argued  that 
torsional  loading  could  become  important  under  some  conditions  in  accounting  for 
potential  asymmetry  associated  with  the  underwater  explosion  bubble  loading  on  external 
hull  appendages  such  as  the  sail  and  planes.  Therefore,  axial  and  torsional  degrees  of 
freedom  will  be  omitted  for  this  analysis,  although  the  entire  3 -dimensional  beam  element 
will  be  used  for  illustration  in  the  following  discussion. 


Figure  3-2.  Two-dimensional,  two-node  beam  finite  element. 
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Figure  3-3.  General  three-dimensional,  two-node  beam  finite  element. 
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In  order  to  simplify  the  assembly  of  the  global  equations  of  motion  for  the  hull,  the 
equations  of  motion  for  each  beam  finite  element  may  be  written  in  the  partitioned  form 
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where  subscripts  1  and  2  refer  to  nodes  1  and  2  of  the  element,  and  each  partition  is  a  6x6 
submatrix  or  6  element  subvector    For  the  genera!  three-dimensional  beam  element,  the 
subvectors  for  nodal  displacement  and  nodal  force  would  be  written  in  the  form: 
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where  subscript  /  refers  to  the  node  ( 1  or  2),  x  is  nodal  translation,  0  is  nodal  rotation,  f  is 
the  (external)  force  acting  on  the  node,  and  m  is  the  (external)  moment  acting  on  the  node 
It  will  be  shown  subsequently  that  formulation  of  submatrices  for  mass  and  stiffness  can  be 
readily  accomplished.  However,  the  definition  of  element  damping  matrices  as  in  equation 
(3-22)  is  generally  not  possible,  thus  a  formulation  of  global  hull  damping  will  be  given  in 
the  next  section  which  will  approximate  the  overall  energy  dissipation  within  the  hull 
structure 

The  stiffness  matrix.  The  derivation  of  stiffness  properties  for  a  finite  element  is 
a  lengthy  process  and  will  not  be  discussed  in  detail  here.   Barltrop  and  Adams77  and 
Hughes78  provide  formulation  for  the  element  stiffness  matrix  of  a  general  three- 
dimensional,  two-noded  beam  finite  element    The  element  stiffness  matrix  can  be  written 
in  a  node-oriented  partitioned  form  consistent  with  equations  (3-22)  and  (3-23).  Each 
partitioned  submatrix  (6x6)  has  terms  ij  which  represent  the  element  elastic  force  in  the 
direction  of  freedom  /  when  a  unit  displacement  is  applied  to  freedom/  with  all  other 
freedoms  held  at  zero  displacement    The  stiffness  submatrices  consistent  with  equation 
(3-22)  can  be  defined  as  shown  in  figure  3-4 

The  mass  matrix.  The  simplest  method  of  assigning  mass  properties  of  the 
ship/beam  model  is  to  assume  that  the  mass  of  each  section  is  lumped  or  concentrated  at 


77Barltrop.  N.D.P..  and  Adams.  A.  J..  Dynamics  of  Fixed  Marine  Structures.  Butterworth-Heinemann. 
UK.  1991.  Appendix  D 
78Hughes.  op.  at.,  pp.  212-217. 
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the  nodal  points  (figure  3-1 )    This  is  called  the  lumped  mass  formulation    In  this  way,  the 
lumped  masses  fill  the  diagonal  terms  of  the  element  mass  matrix  (corresponding  to  the 
individual  nodal  degrees  of  freedom),  while  the  off-diagonal  terms  are  all  zero    The 
matrix  (or  submatrix)  elements  would  therefore  equal  the  mass  of  the  node  for  elements 
corresponding  to  translational  degrees  of  freedom,  or  the  mass  moment  of  inertia  of  the 
node  for  those  elements  corresponding  to  rotational  degrees  of  freedom    The  translation 
of  a  lumped  mass  accounts  for  translational  inertia  effects,  while  the  rotation  of  the 
lumped  mass  accounts  for  rotary  inertia  effects 

The  lumped  mass  formulation  has  some  very  appealing  properties    Firstly,  it  is 
easy  to  associate  a  physical  model  with  the  matrix  (each  degree  of  freedom  at  each  node 
having  its  own  mass).   Secondly,  the  entire  structure  mass  matrix  is  diagonal,  leading  to 
significantly  fewer  computations  and  storage  requirements    The  finite  element  mass 
submatrices  as  defined  in  equation  (3-22),  for  a  lumped  mass  beam  finite  element,  could  be 
written  as  shown  in  figure  3-5. 
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A  is  the  cross-sectional  area  of  the  element,  E  is  Young's  Modulus,  L  is  the  length  of  the 
element,  IZ,IV  are  the  2nd  moments  of  area  of  the  element  cross-section  about  the  z  and 
y  axes,  G  is  the  shear  modulus,  J  is  the  polar  moment  of  inertia  of  the  element  cross- 
section  (about  the  x  axis),  and  A^. ,  Asz  are  the  portions  of  the  cross-sectional  area  over 
which  the  y  and  z-directed  shear  forces  are  assumed  to  act. 
Figure  3-4.  Partitioned  submatrices  for  stiffness  of  a  general  3-D  beam  finite  element. 
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Figure  3-5.  Partitioned  submatrices  for  mass  of  a  general  3-D  beam  finite  element. 

It  is  also  possible  to  take  advantage  of  the  finite  element  formulation  to  develop 
what  is  called  a  consistent  mass  formulation  for  a  beam  element.   In  contrast  to  the 
lumped  mass  matrix,  the  components  of  the  consistent  mass  matrix  cannot  be  associated 
with  a  physical  model.  More  importantly,  after  assemblage,  the  global  mass  matrix  will 
not  be  diagonal,  but  will  include  both  diagonal  and  off-diagonal  terms.  For  this  reason, 
the  consistent  mass  matrix  formulation  requires  considerably  more  computational  effort 
than  the  lumped  mass  formulation.  While  the  consistent  mass  matrix  formulation  is 
generally  superior  in  terms  of  minimizing  errors  in  calculation  of  natural  frequencies79,  the 
errors  associated  with  the  lumped  mass  formulation  significantly  decrease  as  the  number  of 
nodes  used  to  model  the  beam  increases.  Additionally,  when  external  loads  are  being 
applied  directly  to  the  nodes  (as  will  be  discussed  subsequently),  the  equations  of  motion 
behave  better  when  the  inertial  forces  are  likewise  lumped  at  the  nodes.  For  these  reasons, 
the  lumped  mass  matrix  formulation  will  be  used  in  this  submarine  whipping  model. 

3.3.1.2  Hull  damping  model. 

Rayleigh  damping.  As  mentioned  previously,  an  approximation  for  the  global 
damping  matrix  (C)  associated  with  the  hull  structure  can  be  made  by  applying 
appropriate  Rayleigh  coefficients  to  the  global  mass  and  stiffness  matrices  derived  by 
superposition  of  the  element  mass  and  stiffness  matrices: 


79Doyle.  op.  at.,  p.  272. 


C  =  ctM  +  PK  (3-24) 

where  M  and  K  are  the  global  mass  and  stiffness  matrices,  and  coefficients  a  and  J3  are 
the  two  damping  proportionality  constants.   As  discussed  in  section  2.2.3,  coefficients  a 
and  J3  can  be  calculated  if  two  modal  damping  ratios  (£a  and  C,b)  associated  with  two 
specific  modal  frequencies  (co  a  and  co  b)  are  known: 


Q"  -CO 


Wb  "©a 

1/gk     1  /  co 


(3-25) 


The  downside  to  the  Rayleigh  (or  proportional)  damping  matrix  formulation  is  that 
the  contribution  of  the  higher  frequency  modes  of  hull  flexural  response  (modal 
frequencies  much  greater  than  co  b )  are  effectively  eliminated  by  the  high  damping  imposed 
by  the  selection  of  the  particular  coefficients  (see  figure  2-4).  In  the  case  of  elastic  ship 
whipping,  where  only  the  lower  few  modes  of  flexural  response  are  generally  required  to 
adequately  describe  the  motion  of  the  structure,  this  limitation  can  be  considered 
acceptable  as  long  as  the  higher  frequency  used  for  evaluation  of  the  Rayleigh  coefficients 
(co  b)  is  set  among  the  higher  modes  expected  to  contribute  significantly  to  the  whipping 
response.  Thus,  the  problem  becomes  one  of  calculating  the  modal  (natural)  frequencies 
for  the  lowest  modes  of  hull  whipping,  applying  experimentally  determined  (or  reasonably 
approximated)  modal  damping  ratios  for  two  of  the  modes,  and  calculating  the  appropriate 
Rayleigh  coefficients  using  equation  (3-25).  As  discussed  in  section  2.2.3,  the  modal 
(natural)  frequencies  of  a  submarine  structure  can  be  calculated  by  solving  the  generalized 
eigenproblem  for  the  undamped  free  vibration  of  the  dynamic  system 

K<j)n=co;M<j>n  (3-26) 

where  co  n  is  the  natural  frequency  of  vibration  of  the  nth  mode  and  c^  is  the  mode  shape 
vector  of  the  wth  mode.  Data  for  experimentally  determined  (hull)  modal  damping  ratios 
for  the  fundamental  modes  of  submarine  vibrations  is  lacking  in  the  literature.  However, 
some  work  has  been  done  in  the  marine  industry  to  estimate  hysteretic  and  structural 
modal  damping  ratios  for  steel  pile  platforms  for  fundamental-mode  flexural  vibrations  in 
random  seas80,  implying  that  hysteretic  damping  in  thin  steel  cylindrical  structures  is 
typically  less  than  1%  of  critical  (i.e.  C,  <  0.01)  for  the  fundamental  flexural  mode.  There 
is  also  some  older  data  available  in  the  literature  for  modal  damping  ratios  of  steel  circular 
cylinder  ship  masts  based  upon  modal  "bump"  testing81,  suggesting  that  general  steel 


80Cook.  M.F..  and  Vandiver.  J.K..  "Measured  and  Predicted  Dynamic  Response  of  a  Single  Pile  Platform 
to  Random  Wave  Excitation".  Proceedings  of  the  Nth  Offshore  Technology  Conference,  OTC  4285, 
1982.  pp.  637-645. 
8 'Smith,  op.  cit. 
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cylinders  have  (material  and  structural)  modal  damping  ratios  on  the  order  of  2-3  %  of 
critical  for  the  lower  flexural  modes,  although  this  data  seems  less  reliable  because  of 
other  unknown  damping  effects  on  the  mast  not  considered  (air  and  structural  connection 
to  the  ship  hull)    Although  submarines  are  not  precisely  the  same  as  thin  steel  pile 
structures  or  ship  masts,  the  general  hysteretic  damping  mechanisms  are  basically  the  same 
and  thus  it  would  seem  reasonable  to  extend  this  to  data  for  a  similar  cylindrical  steel 
structures  to  a  steel  submarine.  Thus,  it  would  be  reasonable  to  assume  (as  a  first 
approximation)  modal  damping  ratios  (for  the  hull)  on  the  order  of  1%  for  the 
fundamental  flexural  modes  of  a  steel  submarine 

Thus,  once  the  modal  (natural)  frequencies  of  the  hull  structure  have  been 
calculated  by  equations  (3-26),  reasonable  Rayleigh  coefficients  could  be  calculated  by 
equations  (3-25)  assuming  modal  damping  ratios  on  the  order  of  0.01 .  Finally,  the  hull 
damping  matrix  to  be  applied  in  the  finite  element  calculations  can  be  calculated  using  the 
global  mass  and  stiffness  matrices  and  equation  (3-24). 

This  process  will  be  incorporated  into  the  computational  algorithm  and  used  in 
analyses  in  the  next  chapter. 

An  alternative  approach.  "Dry  mode"  modal  superposition.  An  alternative  to 
developing  an  explicit  damping  matrix  (C)  using  Rayleigh  damping  and  solving  the 
complete  equations  of  motion,  is  to  use  the  approximated  hull  modal  damping  ratios 
directly  and  a  modal  superposition  approach.  As  discussed  in  section  3.1,  each  vector 
term  of  the  equations  of  motion  may  be  represented  by  an  independent  formulation  (which 
may  involve  specific  assumptions  of  linearity  for  that  term),  so  long  as  all  the  terms  remain 
coupled  via  a  single  numerical  solution  algorithm.  .Although  forces  associated  with 
external  (fluid)  loading  may  be  nonlinear,  and  the  equations  of  motion  must  therefore  be 
solved  in  the  time  domain,  it  is  still  possible  to  solve  for  the  forces  associated  with  the  hull 
structure  using  a  modal  superposition  technique.  This  type  of  approach  has  been  referred 
to  as  "dry  mode"  modal  superposition82.  In  a  "dry  mode"  modal  superposition  approach, 
all  forces  not  associated  with  the  hull  are  simply  termed  "loading"  forces  and  thus  carried 
on  the  right  hand  side  of  the  equations  of  motion. 

{Mx  +  Cx  +  Kx}hull  =Floading  (3-27) 

Thus,  a  coordinate  transformation  (of  the  generalized  coordinates  of  the  ship  hull)  may  be 
conducted,  so  long  as  the  transformation  of  the  loading  vector  (i.e.  Fmemal  and  Fextemai  as 
discussed  previously)  is  conducted  at  each  time  interval. 


82Bishop.  R.E.D..  and  Pnce.  W.G..  Hydroelasticity  of  Ships.  Cambridge  University  Press.  Cambridge. 
UK.  1979 
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This  formulation  is  consistent  with  the  formulation  of  the  Newmark  time 
integration  method  presented  in  section  3.2  (i.e.  equation  (3-10)).  Because  the  forces 
associated  with  the  hull  are  to  be  modeled  using  assumptions  of  linearity  (as  discussed  in 
the  previous  section),  these  forces  may  be  redefined  using  the  principles  of  modal 
superpostion  outlined  in  section  2.2.3.  The  "dry"  mode  shapes  are  first  found  by  solving 
the  generalized  eigenproblem  for  the  undamped  free  vibration  of  the  hull  structure  using 
equation  (3-26).  The  modal  equations  of  motion  (equations  (2-12))  are  then  written  from 
equations  (3-27)  by  applying  the  modal  transformation,  equations  (2-10)  and  (2-13).  The 
transformation  of  the  modal  force  (i.e.  Qn  in  equations  (2-13))  must  be  performed  at  each 
iteration  at  each  time  step  of  the  direct  time  integration  procedure,  since  the  loading 
forces  (Floading )  are  generally  nonlinear  and  coupled  to  the  generalized  displacements, 

velocities,  and  accelerations  (i.e.  are  functions  of  x,x,x). 

While  this  approach  has  some  appealing  qualities  in  terms  of  representation  of  the 
response  of  a  whipping  ship  in  terms  of  its  fundamental  flexural  mode  shapes,  it  affords 
little  (if  any)  computational  advantage  over  rigorous  computation  using  the  generalized 
coordinates  in  that  two  transformations  must  still  be  made  at  each  iteration  step  of  each 
time  step.  Thus,  the  fundamental  benefit  of  modal  decomposition  is  lost  due  to  the 
requirement  of  maintaining  the  complex  (i.e.  nonlinear)  external  loading  forces  in  terms  of 
the  generalized  coordinate  system.  For  this  reason,  this  approach  will  not  be  pursued  in 
this  analysis  and  all  computation  will  be  carried  out  in  terms  of  generalized  coordinates. 

3.3.2     Internal  forces. 

In  order  to  analyze  the  potential  contribution  of  discrete  internal  structures  and 
sloshing  liquids  on  the  whipping  response  of  a  ship,  it  is  necessary  to  assume  some 
"reasonable"  properties  of  the  internal  structure  and  liquid  storage  tanks  as  might  be 
found  on  a  "real"  ship.    Appropriate  models  to  account  for  the  dynamics  of  the  internal 
structure  and  sloshing  liquid  could  be  made  using  simplified  geometries  and  the  general 
models  discussed  in  chapter  2. 

3.3.2.1  Model  for  discrete  masses. 

The  most  straight-forward  method  of  quantitatively  determining  potential  dynamic 
effects  of  internal  structures  on  the  whipping  response  of  a  ship  is  to  model  the  internal 
structure  with  only  a  few  degrees  of  freedom  (i.e.  modes  of  response),  and  select 
properties  of  the  structure  (i.e.  mass,  stiffness,  damping)  which  are  realistic  and  will  make 
the  internal  structure  respond  near  resonance  (response  modes  not  near  resonance  would 
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contribute  little  to  the  overall  response).  This  amounts  to  finding  a  maximum  effect 
provided  by  an  internal  structure  if  it  were  "tuned"  to  the  resonant  whipping  frequency  (or 
the  exciting  bubble  pulse  frequency).  By  using  this  approach,  a  quantitative  measure  of 
the  damping  effects  of  internal  structures  could  be  obtained  for  various  combinations  of 
mass  sizes  and  locations    While  mass  sizes  should  represent  reasonable  ship  internal 
structures,  the  locations  of  the  masses  may  influence  any  of  the  low  frequency  flexural 
modes  of  the  hull,  and  thus  a  range  of  potential  effects  can  be  predicted.  Figure  3-6 
illustrates  a  notional  internal  structure  connected  to  the  hull  through  some  combination  of 
foundation  and  mounts  with  known  mechanical  properties. 


Figure  3-6.  Notional  internal  structure  connected  to  the  hull 


The  dynamic  vibration  absorber.  For  this  analysis,  a  dynamic  vibration  absorber 
will  be  used  to  model  a  discrete  internal  structure.  The  absorber  is  considered  as  a  single 
discrete  mass,  as  discussed  in  section  2.3. 1,  and  connected  to  the  ship  hull  at  one  of  the 
"nodes"  defined  in  section  3.3.1.   Ideally,  it  would  be  desirable  to  permit  the  absorber  to 
respond  in  all  six  degrees  of  freedom  (three  translational  and  three  rotational).  However, 
in  order  to  reasonably  limit  the  computational  and  modeling  effort,  the  absorber  will  be 
limited  to  respond  only  in  the  vertical  or  horizontal  directions,  as  shown  in  figure  3-7 
This  limitation  simplifies  the  computational  algorithm,  without  unduly  restricting  the  value 
of  the  model.  The  absorber  model  can  be  thought  of  as  being  representative  of  a  resonant 
mode  of  any  discrete  internal  structure. 

Using  the  approach  developed  in  section  2.3. 1  for  a  two  degree  of  freedom 
dynamic  absorber,  the  equations  of  motion  for  the  ship  hull,  as  modified  for  the  forces 
exerted  on  the  hull  by  the  absorber  may  be  written  in  matrix  form: 

[MX  +  CX  +  KX]hun  =  Fextemal  +  Fabsorber  (3-28) 
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where  terms  on  the  left  hand  side  are  defined  by  the  finite  element  model  for  the  ship  hull 
as  described  in  section  3.3.1.  Displacements  and  forces  are  vectors  positioned  at  the  hull 
nodes,  and  hull  mass,  damping,  and  stiffness  are  matrices,  constructed  for  the  finite 
element  formulation.  The  vector  of  forces  exerted  on  the  hull  structure  (i.e.  at  the  nodes) 
by  the  absorber  (Fabsorber )  is  written  (ignoring  damping  associated  with  the  absorber): 

Fab^r  =  {k( x  -  X)}  =  {  - mx}  (3-29) 

where  X  (capital)  refers  to  the  hull  displacement  at  the  node  where  the  absorber  acts,  and 
x  (small  case)  is  the  absorber  displacement  in  the  direction  corresponding  to  X,  as  shown 
in  figure  3-7    Thus,  for  a  finite  element  representation  of  a  ship  hull  with  N  degrees  of 
freedom  (i.e.  N  equations),  one  additional  degree  of  freedom  (i.e.  one  equation)  is 
required  for  each  absorber  used 


K    V  *>         C  A 

^ — s: — k — s — ^ — s — ^ — ^ — k — <: 


Figure  3-7.  Dynamic  vibration  absorber  free  to  respond  in 
vertical  or  horizontal  directions. 


Suppression  of  dynamic  response  of  the  whipping  hull.  As  was  discussed  in 
section  2.3.1,  the  amplitude  of  motion  of  a  structure  could  be  decreased  if  the  properties 
of  a  dynamic  vibration  absorber  were  "tuned"  to  match  the  resonant  response  frequency  of 
the  main  structure.  If  it  is  desired  to  consider  decrease  in  amplitude  of  one  of  the  flexural 
modes  of  hull  whipping  by  application  of  a  dynamic  vibration  absorber  (the  lowest  flexural 
mode  for  example),  then  the  hull  portion  of  the  two  degree  of  freedom  system  depicted  in 
figure  3-7  could  be  represented  by  the  modal  mass,  damping  and  stiffness  for  that 
whipping  mode  (defined  by  equations  (2-13)). 
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In  order  to  quantify  the  potential  dynamic  effects  of  internal  structures  on  the  low- 
frequency  flexural  (whipping)  modes  of  a  submarine  hull,  the  following  process  could  be 
implemented  to  develop  the  equations  of  motion  for  the  forces  exerted  by  the  absorber 
mass  (internal  structure)  on  the  hull: 

1 .  Calculate  the  mode  shapes  and  natural  flexural  frequencies  of  hull  whipping 
(Qn)  by  solving  the  finite  element  eigenproblem  as  discussed  in  section 
3.3.1.  These  represent  the  dry  mode  shapes  and  natural  frequencies  of 
hull  whipping.  For  a  submerged  submarine,  the  mass  matrix  (M)  should  be 
modified  by  adding  approximate  fluid  added  mass  (Ma)  prior  to  calculating 
the  (modal)  natural  frequencies  since  the  added  mass  could  have  a 
significant  effect  on  the  natural  frequencies. 

2.  Calculate  the  modal  mass,  damping,  and  stiffness  for  the  hull  whipping 
mode  of  interest  as  discussed  in  section  2.2.2  (equations  (2-13)). 

3.  Select  a  "reasonable"  absorber  mass  or  mass  of  an  internal  structure  (m) 
(this  could  be  done  for  a  range  of  masses). 

4.  Calculate  the  "optimum"  absorber  stiffness  (k),  which  will  "tune"  the 
absorber  to  the  resonant  frequency  of  the  hull  (or  the  exciting  bubble  pulse 
frequency):   k  =  m-Q~ 

5.  Solve  for  the  dynamic  response  of  the  structure  with  the  internal  absorber 
under  specified  loading  or  initial  conditions  using  equations  (3-28)  and 
(3-29). 

This  process  will  be  incorporated  into  the  computational  algorithm  and  used  in 
analyses  in  the  next  chapter 

3.3.2.2  Model  of  liquid  sloshing  in  a  rectangular  tank. 

Similarly  to  the  case  of  internal  structure,  a  quantitative  measure  of  the  potential 
damping  effects  provided  by  internal  sloshing  liquids  may  be  made  by  considering  the 
forces  exerted  on  the  whipping  hull  by  the  modal  equivalent  mechanical  system  for  the 
complex  internal  sloshing  liquid.  In  effect,  only  those  sloshing  modes  likely  to  be  excited 
by  the  whipping  of  the  hull  need  be  considered  in  the  mechanical  equivalent  system.  The 
sloshing  modes  not  excited  may  be  considered  as  rigid  masses,  and  thus  may  be  considered 
as  part  of  the  main  hull  mass. 

Using  the  relations  developed  in  section  2.3.2,  the  equations  of  motion  for  the  ship 
hull,  as  modified  for  the  forces  exerted  on  the  hull  by  a  sloshing  liquid  may  be  written  in 
matrix  form: 
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[MX  +  CX  +  KX]hui]  =  Fe_al  +  Fhquid  (3-30) 

where  terms  on  the  left  hand  side  are  defined  by  the  finite  element  model  for  the  hull  as 
described  in  section  3.3.1.  Displacements  and  forces  are  vectors  positioned  at  the  hull 
nodes,  and  hull  mass,  damping,  and  stiffness  are  matrices,  constructed  for  the  finite 
element  formulation.  The  vector  of  forces  exerted  on  the  hull  structure  (i.e.  at  the  nodes) 
by  the  sloshing  liquid  (Fh  ld )  can  be  written  (ignoring  damping  in  the  tank): 

F,iquld=-Imax>-m,X  +  ]L[M*n-X)]  (3-31) 

n  =  0  n=l 

where  small  case  letters  refer  to  the  mechanical  equivalent  dynamic  parameters  for  the 
liquid,  and  subscript  /;  refers  to  the  mode  of  sloshing.  Additionally,  the  portion  of  the 
liquid  which  acts  as  a  rigid  body  (m0)  may  be  found  from  equation  (2-46): 

m0=MI-£mn  (3-32) 

n=l 

where  M,  is  the  total  mass  of  liquid  in  the  tank. 

These  relations  may  be  simplified  if  it  is  assumed  that  only  the  fundamental 
sloshing  mode  (n  =  1)  is  excited  near  its  resonance  (the  higher  sloshing  modes  therefore 
may  be  considered  to  act  as  part  of  the  rigid  body  mode,  m  ,).  Equation  (3-31)  and  (3-32) 
may  be  rewritten  and  combined: 

F1,quid  =  -Imnxn  =  -(M1X-m1X  +  m,x1)  =  (m1-M!)X  +  [k1(x,-X)] 

n  =  0 

(3-33) 
Thus,  for  a  finite  element  representation  of  a  ship  hull  with  N  degrees  of  freedom  (i.e.  N 
equations),  one  additional  degree  of  freedom  (i.e.  one  equation)  is  added  for  each  sloshing 
mode  of  interest  and  equations  (3-30)  and  (3-33)  form  a  set  of  N+l  simultaneous 
equations. 

Lateral  sloshing  in  a  rectangular  tank.  Figure  3-8  illustrates  a  rectangular  tank, 
containing  a  liquid  with  a  free  surface,  which  is  oscillating  laterally    The  potential  flow 
solution  for  the  lateral  sloshing  natural  frequencies  and  hydrodynamic  forces  are  given  by 
Vandiver  and  Mitome83  for  the  case  where  horizontal  displacement  of  the  tank  is  assumed 
to  be  harmonic  and  of  the  form: 

X(t)  =  Aeiw'  (3-34) 


83Vandiver  and  Mitome.  op.  at.,  pp.  26-28. 
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where  A  is  the  amplitude  of  the  harmonic  motion  of  the  tank  wall.  The  natural  frequencies 
are  given  for  mode  n: 
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co^  =  g(2n  -  1)  — tanh 
a 


f2n-l 


Tun 


(3-35) 


where  g  is  the  acceleration  due  to  gravity,  a  is  the  width  of  the  rectangular  tank,  and  h  is 
the  mean  depth  of  liquid  in  the  tank    The  modal  masses  are  given  for  each  mode  n: 


m„  =  M, 


8tanhU2n-l)7t 


7U3-(2n-l) 
a 


(3-36) 


As  discussed  in  section  2  3  2,  the  modal  stiffnesses  can  be  computed  from  the  modal  mass 
and  natural  frequency  by  the  expression: 


n       n 


(3-37) 


X(t) 


t 
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Figure  3-8.  Laterally  oscillating  rectangular  tank  with  standing  waves. 

Suppression  of  dynamic  response  of  the  whipping  hull.   Similarly  to  the 
application  of  the  dynamic  vibration  absorber,  the  motion  of  the  whipping  hull  could  be 
analyzed  by  considering  any  onboard  tanks  as  "liquid  slosh"  dampers    The  main 
differences  between  sloshing  liquids  and  dynamic  vibration  absorbers  discussed  in  the 
previous  section  come  about  because  of  the  numerous  sloshing  modes  which  are  not 
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excited  by  the  whipping  (i.e.  the  portion  of  the  liquid  in  the  tank  which  responds  as  a  rigid 
body  with  the  hull,  m0).  Each  sloshing  mode  which  responds  may  be  thought  of  as  an 

individual  vibration  absorber,  but  the  sloshing  modes  which  do  not  respond  might  simply 
be  added  to  the  mass  of  the  hull  structure. 

Thus,  a  quantitative  feel  for  the  effects  of  sloshing  liquids  might  be  obtained  by 
considering  the  effects  of  internal  structures  or  absorbers,  but  also  considering  the  portion 
of  the  liquids  which  are  not  excited,  but  respond  as  rigid  masses.  A  more  thorough 
discussion  of  this  will  be  presented  in  the  next  chapter. 

3.3.3     External  (fluid)  forces. 

3.3.3.1  Representation  of  the  hull  surface. 

As  discussed  in  section  2.3.1,  slender  bodies  are  often  represented 
hydrodynamically  using  a  "strip  theory"  approach  since  flow  around  each  strip  of  the  hull 
tends  to  be  represented  well  by  considering  the  flow  to  be  two-dimensional  (particularly  at 
higher  frequencies).  Slender  submarine  hulls,  whose  length-to-diameter  ratios  typically 
range  from  eight  to  twelve,  would  therefore  be  represented  well  using  the  hydrodynamic 
"strip  theory"  approach. 

For  purposes  of  this  analysis,  a  submarine  hull  will  be  represented 
hydrodynamically  using  a  number  of  surface  elements  ("strips")  as  shown  in  figure 
3-9.  The  centers  of  pressure  of  the  surface  elements  will  be  considered  to  be  located 
coincident  with  the  beam  finite  element  nodes  defined  in  section  3.3.1.1.  Thus,  the  net 
forces  exerted  on  each  of  the  hull  sections  by  the  surrounding  fluid  could  be  considered  to 
act  at  the  nodes  of  the  finite  element  model.  The  properties  of  each  surface  element  (i.e. 
principle  dimensions)  would  be  considered  constant  for  the  length  of  that  surface  element. 

Because  two-dimensional  flow  is  assumed  around  each  surface  element,  the  forces 
on  each  surface  element  (and  therefore  on  each  finite  element  node)  will  generally  be 
resolved  into  two  Cartesian  coordinate  components  corresponding  to  the  horizontal  and 
vertical  directions  (normal  to  the  ship  axis).  Force  components  on  each  surface  element 
will  be  computed  and  applied  to  the  beam  finite  element  node.  In  the  case  where  a 
particular  surface  element  has  appendages,  additional  contributions  to  the  horizontal  and 
vertical  forces  will  be  computed  for  the  given  appendage  geometry. 
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Surface   element   ('strip') 


Figure  3-9.   Hydrodynamic  representation  of  the  submarine  hull  by 
surface  elements  ("strips") 


3.3.3.2  Morison  relative  velocity  model. 

As  discussed  in  section  2.3.2,  when  two  dimensional,  generally  oscillating  flow  can 
be  assumed,  and  when  both  fluid  inertia  and  drag  forces  are  important  (and  wave  radiation 
may  be  ignored),  then  a  Morison's  equation  may  be  used  to  predict  hydrodynamic  forces 
on  a  body    It  will  be  assumed  for  this  analysis,  that  the  whipping  submarine  is  sufficiently 
below  the  free  surface  such  that  wave  radiation  may  be  ignored    Under  this  assumption, 
the  hydrodynamic  forces  on  each  of  the  surface  elements  may  be  predicted  using  the 
Morison  relative  velocity  formulation  given  by  equation  (2-62)  and  corrected  by 
multiplying  by  the  length  of  the  element: 


hvdro 


(  i  i       A 

pACm(u- x)  +  pAx  +  — pDCd  (u-x)  u-x    L 


(3-38) 


where  Fhvdril  is  the  total  hydrodynamic  force  on  the  surface  element,  p  is  the  fluid  density, 

A  is  the  cross-sectional  area  of  the  surface  element  (normal  to  the  direction  of  flow),  D  is 
the  characteristic  diameter  of  the  surface  element  (normal  to  the  direction  of  flow),  L  is 
the  length  of  the  surface  element,  u  is  the  fluid  velocity  at  the  location  of  the  center  of 
pressure  of  the  surface  element  (i.e  the  "free  field"  fluid  velocity  caused  by  the  pulsating 
explosion  bubble  and  considered  relatively  unaffected  by  the  presence  of  the  hull),  and  x  is 
the  displacement  of  the  center  of  pressure  of  the  surface  element  (i.e.  the  finite  element 
nodal  displacement)    Cm  and  Cd  are  the  relative  motion  fluid  inertia  and  drag 
hydrodynamic  coefficients  (respectively)  and  are  defined  by  equation  (2-63) 

Determination  of  hydrodynamic  coefficients.  For  the  case  of  the  whipping  of  a 
submerged  submarine,  particular  care  must  be  taken  to  select  hydrodynamic  coefficients 
for  analysis  which  properly  represent  both  the  geometry  and  the  conditions  of  the  relative 


fluid  flow  around  the  surface  element.  Generally,  it  would  be  expected  that  separation  of 
flow  around  the  hull  would  occur  during  submarine  whipping  motions,  particularly  around 
those  sections  with  small  characteristic  dimensions  (or  containing  components  with  small 
characteristic  dimensions).  Determination  of  the  extent  to  which  separation  occurs  and 
the  resulting  forces  generally  requires  a  complex  evaluation  of  the  separation  and  vortex 
shedding  mechanisms    In  the  "discrete-point"  vortex  method84,  for  example,  forces  on  an 
isolated  sharp  edge  may  be  predicted  my  modeling  the  separating  shear  layers  which  feed 
the  growing  vortices.  The  hydrodynamic  coefficients  are  calculated  iterativelv  over  a 
number  of  cycles  by  a  Fourier  averaging  process    This  type  of  evaluation,  while  perhaps 
required  for  completely  rigorous  solution  of  the  hydrodynamic  forces,  is  beyond  the  scope 
of  this  analysis  and  will  not  be  considered  here. 

An  initial  perspective  of  the  flow  fields  which  could  be  expected  around  a 
whipping  submarine,  and  thus  of  the  criteria  necessary  for  initial  selection  of 
hydrodynamic  coefficients  for  analysis,  could  be  attained  by  considering  past  experience 
with  ship  whipping,  and  considering  typical  submarine  geometries.  As  pointed  out  by 
Hicks85,  the  fundamental  (i.e.  lowest)  natural  whipping  period  of  most  ships  is  typically 
between  0.3  and  1 .0  seconds  (fairly  well  matched  with  typical  bubble  pulse  periods).  For 
typical  characteristic  dimensions  (D)  of  submarine  hull  sections  (ranges  might  be  between 
0. 1  m  for  small  control  surfaces,  to  15  meters  for  an  entire  hull  section),  characteristic 
frequencies  (3)  would  range  from  1  x  104  to  7.5  x  108  (higher  whipping  modes  would 
yield  even  higher  3)    These  characteristic  frequencies  correspond  to  very  large  Reynolds' 
numbers  (Rn)  over  a  large  range  of  Keulegan-Carpenter  numbers  (KC). 

Extensive  experimental  studies  have  been  conducted  in  the  marine  industry  in 
order  to  obtain  values  for  hydrodynamic  coefficients  Cm,  Cd,  CA,  and  Cj,  particularly  for 
smooth  horizontal  cylinders,  as  discussed  in  section  2.3.2.  Even  with  the  extensive 
studies,  a  great  deal  of  experimental  "scatter"  exists  in  the  data  base,  due  to  such  factors 
as  free  surface  effects,  surface  roughness,  and  even  the  random  nature  and  directionality  of 
the  ocean  environment  (for  field  tests).  For  these  reasons,  a  great  deal  of  care  must  be 
taken  prior  to  applying  coefficient  values  from  the  literature  to  a  design  or  analysis  case. 
Additionally,  most  experimental  studies  have  been  carried  out  at  relatively  low  frequencies 
(and  therefore  low  Reynolds'  numbers),  and  thus  there  is  limited  data  which  could  be 
directly  applied  to  submarine  whipping  without  some  extrapolation. 


84Graham.  J.M.R..  "The  forces  on  sharp-edged  cylinders  in  oscillatory  flow  at  low  Keulegan-Carpenter 
numbers" .  Journal  of  Fluid  Mechanics.  97.  1980.  pp.  331-346. 
85 Hicks,  op.  cit..  p.  393. 
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Extensive  laboratory  research  conducted  by  Sarpkaya86  showed  that,  for  smooth 
circular  cylinders  at  high  Reynolds'  numbers  (and  therefore  high  (3),  the  value  of  Cd 
approaches  0.65  and  Cm  approaches  1.8  (relatively  constant  over  a  wide  range  of  KC). 
These  coefficients  should  provide  good  first  approximations  for  inertia  and  drag 
coefficients  of  the  circular  hull  sections  (surface  elements)  for  high-frequency  whipping. 

Generally,  there  may  be  several  sharp-edged  portions  of  a  submarine  surface 
(control  surfaces  and  sail)  which  may  induce  significant  vortex  shedding  and  therefore 
significant  viscous  drag  on  particular  hull  surface  elements    Some  accounting  must 
therefore  be  made  for  the  additional  forces  exerted  on  these  surface  elements    Test  data 
on  sharp-edged  cylinders  in  high  Reynolds'  number  oscillatory  flow  is  generally  lacking  in 
the  literature.   Some  data  is  presented  by  Bearman,  et.  al.87  for  flat  plates  in  oscillatory 
flow  at  Keulegan-Carpenter  numbers  between  1  and  10  (relatively  low)  which  might  be 
extrapolated  to  predict  that  Cd  would  fall  in  the  range  of  4-6  and  Cm  would  fall  in  the 
range  of  around  1.0  for  higher  Reynolds'  number  (higher  frequency)  flows    While  these 
values  are  probably  not  experimentally  sound,  they  should  at  least  provide  reasonable  first 
order  approximations  of  additional  viscous  forces  exerted  on  certain  hull  surface  elements 
by  sharp-edged  control  surfaces  or  sail. 

The  modeling  of  hydrodynamic  forces  on  a  whipping  submarine  using  equation 
(3-38)  will  be  included  as  part  of  the  computational  algorithm  developed  in  section  3.4. 
One  analysis  application  will  be  made  in  chapter  4  where  whipping  response  of  the 
submarine  is  due  only  to  an  initial  displacements  (i.e.  free  whipping),  and  therefore  no 
"free  field"  fluid  velocity  or  acceleration  terms  (i.e  u  and  li )  will  be  included,  and  the 
hydrodynamic  forces  will  be  due  only  to  the  motion  of  the  hull  in  otherwise  still  water. 
Additional  analyses  will  include  cases  where  the  submarine  is  excited  by  a  pulsating 
explosion  bubble,  and  "free  field"  fluid  motion  will  be  included  in  the  hydrodynamic 
model. 

3.3.3.3  Model  of  fluid  loading  from  a  pulsating  explosion  bubble. 

In  analysis  cases  where  it  is  desirable  to  couple  ship  whipping  with  the  loading 
imposed  by  the  pulsation  of  an  explosion  bubble,  it  is  necessary  to  define  a  time-varying 
formulation  for  the  "free  field"  fluid  velocity  and  acceleration  associated  with  the  pulsating 


86Sarpkaya.  T..  "In-line  and  transverse  forces  on  cylinders  in  oscillaung  flow  at  high  Reynolds'  numbers". 
Proceedings  of  the  Eight  Offshore  Technology  Conference.  OTC  2533.  1976.  pp.  95-108. 
87Bearman.  P.W..  Downie.  M.J..  Graham.  J.M.R..  and  Obasaju.  ED.  "Forces  on  cylinders  in  viscous 
oscillator*  flow  at  low  Keulegan-Carpenter  numbers" .Journal  of Fluid  Mechanics.  154.  1985.  pp.  344- 
345. 
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bubble  (i.e.  u  and  u )    Because  derivation  of  equations  for  bubble  dynamics  is  a  lengthy 
process  and  not  a  goal  of  this  analysis,  formulation  available  in  the  literature  will  be  only 
briefly  summarized    A  particularly  useful  formulation  for  spherically  symmetric,  pulsating 
and  migrating  explosion  bubbles  is  presented  by  Wilkerson88,  which  is  based  upon  the 
published  theories  of  Herring89,  Friedman90,  Taylor91,  Cole92,  and  Hicks93    This 
formulation  includes  accounting  of  effects  due  to  the  presence  of  a  free  surface  and  drag 
forces  on  the  migration  of  the  bubble. 

The  free  field  equations  of  motion  for  a  spherically  symmetric,  pulsating  and 
migrating  explosion  bubble  are  described  adequately  using  inviscid,  irrotational  fluid  flow 
theory  (i.e.  potential  flow).   Figure  3-10  illustrates  an  appropriate  coordinate  system  which 
encompasses  the  explosion  bubble,  its  "image  bubble"  (to  account  for  the  presence  of  the 
free  surface),  and  an  arbitrary  point  in  the  fluid.  The  fluid  velocity  potential  function  (§) 
resulting  from  a  pulsating  and  migrating  explosion  bubble  may  be  written  as  the 
superposition  of  potentials  resulting  from  a  simple  source  and  its  image  (representing 
pulsation  of  the  spherical  bubble)  and  potentials  resulting  from  a  dipole  and  its  image 
(representing  migration  of  the  spherical  bubble) 

e.      e^  e,      e^ 

(j)  =  — L-r-^cos91 — -  +  ^-cos0,  (3-39) 

where  e,  and  e,  are  the  time-dependent  strengths  of  the  sources  and  dipoles 
(respectively),  r,  ,r,  are  the  radial  distances  from  the  center  of  the  bubbles  to  the  arbitrary 
point  in  the  fluid,  and  0,  ,0,  are  the  angles  from  the  vertical  (as  shown  in  figure  3-10). 
The  first  term  on  the  right  hand  side  is  the  potential  due  to  the  source  representing  the 
pulsation  of  the  explosion  bubble.  The  second  term  is  the  potential  due  to  the  dipole 
representing  the  migration  of  the  explosion  bubble    The  last  two  terms  represent  the 
potentials  due  to  the  effects  of  the  "image  bubble"  which  accounts  for  the  effect  of  the  free 
surface  (i.e.  satisfies  the  boundary  condition  at  the  free  surface). 


88  Wilkerson.  S.A..  "Elasuc  Whipping  Response  of  Ships  to  an  Underwater  Explosion  Loading".  Master  of 

Science  in  Engineering  Thesis.  George  Washington  University.  1985. 

89Herring.  C,  "Theory  of  the  Pulsations  of  the  Gas  Bubble  Produced  by  an  Underwater  Explosion". 

Underwater  Explosions  Research,  I  bl.  II  -  The  Gas  Globe.  Office  of  Naval  Research.  1950. 

90Friedman.  B..  "Theory  of  the  Underwater  Explosion  Bubbles".  Underwater  Explosions  Research,  Vol.  II 

-  The  Gas  Globe.  Office  of  Naval  Research.  1950. 

91Taylor.  G.I..  "Vertical  Mouon  of  a  Spherical  Bubble  and  the  Pressure  Surrounding  It".  Underwater 

Explosions  Research,  Vol.  II  -  The  Gas  Globe.  Office  of  Naval  Research.  1950. 

92Cole.  R.H..  Underwater  Explosions.  Princeton  University  Press.  Princeton.  NJ.  1948. 

93Hicks.  A.N..  "The  Theory  of  Explosion-Induced  Ship  Whipping".  Naval  Construction  Research 

Establishment.  XCRE  Report  R579.  1972. 
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Image"    bubble 


Uy  .Liy  (vert  fluid   vel.acc) 


u„  .(J,  (horiz   fluid    vel.acc) 


Explosion    bubble 


Figure  3-10    Coordinate  system  for  explosion  bubble  dynamics. 
Explosion  and  "image"  bubble  representation. 

The  time-dependent  source  and  dipole  strengths  are  those  required  to  maintain 
continuity  of  the  normal  fluid  particle  velocity  at  the  bubble  surface  (i.e.  the  normal  fluid 
particle  velocity  must  equal  the  rate  of  change  of  the  bubble  radius,  a  ,  corrected  for  the 
migration  velocity).  The  required  source  and  dipole  strengths  are: 


e,  =  a~a 


3  / 


e. 


a"a 
4^ 


(3-40) 


where  vm  is  the  upward  velocity  of  the  (migrating)  bubble  and  d  is  the  depth  of  the 
bubble.  Using  conservation  of  energy,  the  differential  equations  which  govern  the 
pulsation  and  migration  of  the  bubble  (i.e.  solve  for  the  bubble  radius  and  location)  may  be 
derived94.  The  differential  equations  (in  non-dimensionalized  form)  which  describe  bubble 


94Wilkerson.  op.  cit..  pp.  7-24. 
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pulsation  and  migration  are  presented  as  part  of  figure  3-11,  along  with  other  defined 
constants  and  non-dimensionalized  variables  used  in  the  solution  of  the  bubble  dynamics. 
The  first  order  differential  equations  which  govern  the  pulsation  and  migration  of 
the  bubble  may  be  integrated  using  a  time-stepping  integration  procedure  such  as  a  Runge- 
Kutta  method.  It  is  recommended95  that  a  variable  time-stepping  procedure  be  used 
because  of  the  rapid  changes  in  radius  during  the  phases  of  bubble  minima.  With  the 
calculated  bubble  radii  and  positions  (with  time),  the  velocity  potential  may  be  calculated 
using  equations  (3-39)  and  (3-40),  and  the  "free  field"  fluid  velocity  for  any  arbitrary  point 
in  the  fluid  may  be  calculated  by  differentiating  the  velocity  potential  (equation  (2-47)). 
Wilkerson96  derived  expressions  for  fluid  velocities  and  accelerations  (horizontal  and 
vertical  components)  at  an  arbitrary  point  in  the  fluid  as  functions  of  position  and  source 
and  dipole  strengths.   The  resulting  expressions  (given  in  Cartesian  coordinates  of  figure 
3-10)  are  given  in  figure  3-12. 


95Hicks.  op.  at.,  p.  159. 
96Wilkerson.  op.  at.,  pp.  25-29. 
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Differential  equations  which  govern  bubble  pulsation  and  migration: 
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Initial  Conditions: 
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Non-dimensionalized  defined  variables:      a  =  a/L       (non-dimensional  bubble  radius) 

;  =  yb  /  L       (non-dimensional  bubble  depth) 
i  =  t  /  T      (non-dimensional  time) 


where:  L  *  13.6(W/D0)     for  TNT 


k«(0.0552)D225    for  TNT 


T  =  2.94 


W 


1  3 


TV  ° 


for  TNT 


Y  *  125   for  TNT 


C    =  ?  75 


D0=d  +  33         d  is  initial  depth  of  charge       W  is  charge  weight  of  TNT 
Control  variables:         e  (0  for  no n- migrating;  1  for  migrating) 

J3  (0  for  no  free-surface  effect;  1  for  free-surface  effect) 

Figure  3-11.  Non-dimensionalized  variables  and  differential  equations  which  describe 
the  pulsation  and  migration  of  an  underwater  explosion  bubble. 
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where  source  and  dipole  strengths  may  be  defined  (in  non-dimensionalized  form): 


e,  =  Ua'a  IT 


— r(a:d  +  2aa:) 


e.  = 


2T 


L4 


a3X.  +  3ct/jj 


T-  v  2T- 

and:         vm  =  -£L/T  =  -AL/T      (the  upward  velocity  of  the  bubble) 


Figure  3-12.  Fluid  velocities  and  accelerations  (in  relative  Cartesian  coordinates)  for 
given  source  and  dipole  strengths  and  upward  bubble  velocity. 
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3.4        Formulation  of  the  computational  algorithm. 

In  order  to  integrate  each  of  the  desired  force  components  presented  in  the 
previous  sections,  a  comprehensive  computational  algorithm  must  be  developed.  The 
required  algorithm  must  be  capable  of  performing  the  following: 

1.  Compute  the  time  histories  of  the  explosion  bubble  radius  and  depth  by 
integration  of  the  differential  equations  governing  bubble  dynamics  using  a  Runge- 
Kutta  time  integration  procedure. 

2.  Compute  the  "free-field"  fluid  velocities  and  accelerations  due  to  the 
pulsating  and  migrating  explosion  bubble.  These  velocities  and  accelerations  must 
be  calculated  at  the  locations  of  each  of  the  finite  element  nodes  used  to  model  the 
dynamics  of  the  hull  structure  at  each  iteration  at  each  time  step,  and  must  be 
resolved  into  the  vertical  and  horizontal  directions  defined  by  the  ship  coordinate 
system  (a  change  of  coordinates  from  the  bubble  coordinate  system  to  the  ship 
coordinate  system  is  required). 

3.  Compute  the  external  and  internal  loading  on  the  ship  hull  using  the 
appropriate  models  developed  in  the  previous  sections. 

4.  Iteratively  determine  the  hull  structure  displacements,  velocities,  and 
accelerations  using  the  Newmark  method  (with  Newton-Raphson  iterations 
conducted  at  each  time  step)  as  described  in  section  3.2.  Also,  iteratively 
determine  displacements,  velocities,  and  accelerations  of  internal  masses  as 
required. 

A  computational  algorithm  satisfying  the  above  requirements  may  be  generated 
using  a  suitable  computer  language,  separating  the  algorithm  into  sections  or  subroutines 
appropriate  for  computation  using  the  chosen  language    For  this  application,  the 
algorithm  will  be  implemented  on  a  personal  computer  using  subroutines  written  for  the 
MATLAB®  computer  program.   Appendix  A  presents  MATLAB®  subroutines  written 
for  this  analysis,  each  of  which  addresses  specific  portions  of  the  overall  algorithm.  The 
subroutines  are  presented  as  follows: 

a.  Subroutine  IN.M  -  Serves  as  an  input  file  for  other  subroutines.  Data  is  loaded 
from  formatted  batch-type  files.  Arrays  and  constants  used  in  other  subroutines  are 
generated. 

b.  Subroutine  HULL.M  -  Calculates  global  mass  and  stiffness  matrices  for  the  hull 
structure,  solves  the  generalized  eigenproblem  to  obtain  the  "dry  mode"  natural 
frequencies  and  mode  shapes,  calculates  an  appropriate  viscous  hull  damping  matrix  for 
the  hull  using  a  Rayleigh  damping  approach,  and  calculates  vectors  for  modal  mass, 
damping,  and  stiffness. 
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c.   Subroutine  WETMODE.M  -  Calculates  an  approximate  fluid  added  mass 
matrix  for  the  hull  using  a  Morison  equation  approach  (i.e.  strip  theory),  adds  the  fluid 
added  mass  matrix  to  the  dry  hull  mass  matrix  and  solves  the  generalized  eigenproblem  to 
obtain  the  "wet  mode"  natural  frequencies  and  mode  shapes,  and  calculates  vectors  for 
"wet  mode"  modal  mass,  damping  and  stiffness. 

d    Subroutine  FREE  WHIP  M  -  Calculates  nodal  point  response  (displacements, 
velocities,  and  accelerations)  during  free  vibration  with  initial  displacements  given  by 
scaled  mode  shape  vectors  for  specified  modes    Response  is  calculated  using  the 
Newmark  time  integration  method  with  iteration  conducted  at  each  time  step  using  a 
Newton-Raphson  method.  This  subroutine  is  designed  primarily  to  compare  damping 
associated  with  hull  mechanisms  (i.e.  material  hysteretic  and  structural)  and  damping 
associated  with  external  hydrodynamic  forces. 

e.  Subroutine  BUBBLE  M  -  Computes  and  saves  the  time  histories  of  the 
explosion  bubble  radius  and  bubble  depth  by  integration  of  the  differential  equations 
governing  bubble  dynamics    Time  histories  are  computed  using  a  Runge-Kutta  integration 
method  with  a  variable  time  step.  The  differential  equations  governing  bubble  radius  are 
called  from  function  routines  RAD  FN  M  and  RFUNC.M    The  differential  equations 
governing  bubble  depth  are  called  from  function  routines  DPTHFN.M  and  DFUNC.M. 

f.  Subroutine  BUBWHIP.M  -  Calculates  the  nodal  point  response  (displacements, 
velocities,  and  accelerations)  during  forced  vibration  due  to  a  pulsating  explosion  bubble. 
Response  is  calculated  using  the  Newmark  time  integration  method  with  iteration 
conducted  at  each  time  step  using  a  Newton-Raphson  method.  Additionally,  local  axial 
strains  are  calculated  at  specified  hull  locations  using  simple  beam  theory  (used  for 
comparison  of  the  computational  algorithm  with  model  test  data).  The  free-field  fluid 
velocities  and  accelerations  due  to  the  pulsating  explosion  bubble  are  calculated  at  each 
iteration  when  function  routine  FLUID  M  is  called  from  within  BUBWTTIP  M    This 
subroutine  is  designed  primarily  to  provide  a  comprehensive  computation  of  explosion- 
induced  whipping  effects    Additionally,  it  may  be  used  to  quantify  damping  associated 
with  internal  structures  and  to  provide  for  comparison  of  the  computational  algorithm  with 
model  test  data. 


4.  ANALYSIS  OF  WHIPPING  OF  A  SUBMERGED  SUBMARINE 

4.1        Introduction. 

In  this  chapter,  the  prediction  of  the  whipping  response  of  a  submerged  submarine 
is  discussed,  including  an  analysis  of  specific  damping  mechanisms    The  computational 
algorithm  developed  in  the  previous  chapter  is  utilized  to  predict  the  response  of  a 
notional  7400  long  ton  fast  attack  submarine  (SSN).  Additionally,  the  computational 
algorithm  is  used  to  predict  the  response  of  the  U.S.  Navy  test  platform  "Red  Snapper", 
subjected  to  specified  explosive  charge  weight/standoff  combinations  for  comparison  with 
model  test  data.  The  latter  analysis  is  presented  in  non-specific  response  units  due  to 
security  concerns  with  the  actual  test  data  used  for  comparison 

General  characteristics  of  the  notional  submarine  (SSN)  are  given  in  Table  4-1. 
From  these  general  characteristics,  specific  characteristics  (beam  element  data,  hull  section 
data,  etc.)  are  formulated  and  provided  as  input  via  subroutine  IN.M.  The  "dry"  mode 
natural  frequencies  and  mode  shapes  are  calculated  using  subroutine  HULL.M,  and  the  . 
"wet"  mode  natural  frequencies  and  mode  shapes  are  calculated  using  subroutine 
WETMODE.M  (i.e.  the  effects  of  fluid  added  mass  are  included)    Table  4-2  gives  both 
"dry  mode"  and  "wet  mode"  natural  frequencies  computed  for  the  first  12  displacement 
modes.   It  should  be  noted  that  this  represents  the  first  6  displacement  modes  in  each 
plane  (i.e.  horizontal  and  vertical),  of  which  the  first  2  are  considered  rigid  body  modes 
(corresponding  to  zero  natural  frequencies),  and  the  remaining  4  modes  are  considered  as 
the  fundamental flexural  modes.  The  horizontal  "wet"  mode  shapes  are  plotted  (as 
normalized  nodal  displacements)  and  shown  in  figure  4-1  (rigid  body  mode  shapes)  and 
figure  4-2  (fundamental  flexural  mode  shapes). 


Table  4-1 .  General  characteristics  of  the  notional  SSN. 


Length  (ft) 

350 

Displacement/weight  (long  tons) 

7440 

Maximum  beam/diameter  (ft) 

35 

Hull  material 

HY-80  steel 

Number  of  hull  sections/beam  elements 

20/19 

Appendages 

Sail,  bow  planes,  stern  planes,  rudder 

Table  4-2.  Ri 

gid  body  and  fundamental  flexural  "dry"  mode  and  "wet"  mode  natural 
frequencies  for  notional  SSN. 

Mode 

"Dry  mode"  natural  frequency,  co 
(hz,  rad/sec) 

"Wet  mode"  natural  frequency,  co 
(hz,  rad/sec) 

1  Rigid  body 

0 

0 

2  Rigid  body 

0 

0 

3  Rigid  body 

0 

0 

4  Rigid  body 

0 

0 

5  Flexural  (V) 

1.984,    12  466 

1.370,   8.608 

6Flexural  (H) 

1  984,    12.466 

1.373,  8627 

7  Flexural  (V) 

4.498,  28.262 

3  060,   19227 

8  Flexural  (H) 

4498,  28  262 

3.074,   19.315 

9  Flexural  (V) 

7.054,  44.322 

4.827,  30.329 

lOFlexural(H) 

7.054,  44.322 

4.863,  30.555 

11  Flexural  (V) 

9802,  61.588 

6907,  43.398 

12Flexural(H) 

9802,  61.588 

6922,  43.492 

Wet  mode  shape  of  rigid  body  mode  in  horizontal  (y=0)  plane 


2  4  6  8  10         12         14         16         18 

Wet  mode  shape  of  rigid  body  mode  in  horizontal  (y=0)  plane 


2  4  6  8  10         12         14         16         18 

Figure  4-1 .  Horizontal  rigid  body  "wet"  mode  shapes  for  notional  SSN. 
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Wet  mode  shape  of  flexural  mode  1  in  horizontal  (y=0)  plane 


2  4  6  8  10         12         14        16        18 


Wet  mode  shape  of  flexural  mode  2  in  horizontal  (y=0)  plane 


2  4  6  8  10         12         14        16        18 


Wet  mode  shape  of  flexural  mode  3  in  horizontal  (y=0)  plane 


2  4  6  8  10         12         14        16        18 


Wet  mode  shape  of  flexural  mode  4  in  horizontal  (y=0)  plane 


2  4  6  8  10         12         14        16        18 


Figure  4-2.  Horizontal  fundamental  flexural  "wet"  mode  shapes  for  notional  SSN. 
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4.2        Prediction  of  "dry"  whipping  response.  Hull  damping. 

The  "dry"  whipping  response  of  the  SSN  is  computed  for  given  initial  conditions 
using  subroutine  FREE  WHIP  M  (with  fluid  density  set  to  zero  in  Morison's  equation). 
For  free  "dry"  whipping  (as  well  as  free  "wet"  whipping  discussed  subsequently),  initial 
nodal  point  displacements  are  specified  as  scaled  horizontal  mode  shape  vectors,  and  time 
histories  of  nodal  point  displacements,  velocities,  and  accelerations  are  computed  based 
upon  these  initial  displacements    Thus,  initial  loading  en  the  hull  structure  is  due  onlv  to 
internal  forces  due  to  structural  stiffness. 

For  this  analysis,  initial  displacements  are  set  to  the  first  and  second  horizontal 
fundamental  flexural  mode  shapes  (scaled  by  2%  of  the  total  ship  length).  Appendix  B 
presents  plots  of  horizontal  displacements,  velocities,  and  accelerations  for  nodes  7,  10, 
and  20  for  initial  displacements  (a)  scaled  from  the  1st  horizontal  fundamental  flexural 
mode  shape  (i.e.  mode  6  in  Table  4-2),  and  (b)  scaled  from  the  2nd  horizontal  fundamental 
flexural  mode  shape  (i.e.  mode  8  in  Table  4-2). 

As  discussed  in  chapter  3,  it  has  been  assumed  that  the  inherent  hull  damping  (i.e. 
due  to  material  hysteresis  and  structural  damping)  would  reasonably  follow  experimental 
data  for  steel  cylindrical  marine  structures  and  ship  masts.  Thus,  it  was  assumed  that 
modal  damping  ratios  on  the  order  of  1%  for  hull  damping  would  be  sufficient  and 
Rayleigh  coefficients  could  be  calculated  on  this  basis.  Thus,  reductions  in  free  response 
amplitude  with  time  (damping)  seen  in  graphs  presented  in  Appendix  B  are  merely 
manifestation  of  this  initial  assumption,  and  no  further  conclusions  may  be  drawn 
concerning  the  validity  of  this  initial  assumption. 

However,  an  approximate  graphical  check  could  be  performed  to  verify'  that 
damping  in  the  case  of  "dry"  whipping  is  indeed  represented  computationally  by  the 
assumed  1%  modal  damping.  Indeed,  measurement  of  amplitude  log  decrements  (recall 
equation  (1-1))  and  calculation  of  modal  damping  ratios  (recall  equation  (2-40))  for  both 
mode  1  horizontal  flexure  (represented  by  node  10  displacement)  and  mode  2  horizontal 
flexure  (represented  by  node  7  displacement),  show  that  graphically  determined  modal 
damping  ratios  are  approximately  1%. 
4.3        Prediction  of  "wet"  whipping  response.  Hydrodynamic  damping. 

The  "wet"  whipping  response  of  the  SSN  is  computed  for  the  same  given  initial 
conditions  as  the  "dry"  mode  analysis,  also  using  subroutine  FREEWHIP  M  (with 
inclusion  of  hydrodynamic  forces  by  Morison's  equation).  Appendix  C  presents  plots  of 
horizontal  displacements,  velocities,  and  accelerations  for  nodes  7,  10,  and  20  for  initial 
displacements  (a)  scaled  from  the  1st  horizontal  fundamental  flexural  mode  shape,  and  (b) 
scaled  from  the  2nd  horizontal  fundamental  flexural  mode  shape. 
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Figure  4-3  shows  comparisons  between  nodal  point  displacements  for  "wet"  and 
"dry"  cases  for  initial  displacements  scaled  from  mode  1  and  mode  2  flexural  mode  shapes 
Even  though  modal  damping  cannot  be  directly  applied  to  the  nonlinear  hydrodynamic 
damping  forces,  a  calculated  "equivalent"  modal  damping  might  still  serve  as  a  measure  of 
hydrodynamic  damping  in  comparison  to  the  "dry"  case.  As  was  done  for  "dry"  whipping, 
a  graphical  approximation  of  "equivalent"  modal  damping  ratios  can  be  made  for  mode  1 
horizontal  flexure  and  mode  2  horizontal  flexure  (represented  by  node  10  and  node  7 
displacements,  respectively)    As  expected,  the  inclusion  of  external  hydrodynamic  forces 
(represented  by  Morison's  equation)  provide  for  increased  modal  damping  ratios  for  both 
flexural  modes.  Both  mode  1  and  mode  2  flexure  are  characterized  by  "equivalent"  modal 
damping  ratios  of  approximately  2.2%. 

Horizontal  displacement,  node  10 

Dry      Wet 


0 


0.5  1  1.5  2 

Time  (sec) 
Horizontal  displacement,  node  7 


2.5 


1.5 
Time  (sec) 

Figure  4-3.  Comparison  between  nodal  point  displacements  for  "dry"  and  "wet"  hull 

flexure.  Top:  mode  1  flexure,  represented  by  node  10  displacement 

Bottom:  mode  2  flexure,  represented  by  node  7  displacement. 
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4.4        Prediction  of  whipping  response  due  to  explosion  bubble  pulse  loading. 

The  previous  sections  presented  analysis  based  upon  an  assumption  of  initial 
displacements  set  to  flexural  mode  shapes  (i.e.  free  whipping),  aimed  only  at  providing 
some  quantitative  measure  of  hydrodynamic  damping.  As  discussed  in  chapter  1,  the  main 
concern  with  submarine  whipping  involves  the  quasi-periodic  bubble  pulse  loading 
resulting  from  near  to  intermediate  stand-off  explosions    It  is  this  quasi-periodic  bubble 
pulse  loading  which  can  reinforce  the  free  vibration  and  magnify  the  whipping  response, 
particularly  when  the  period  of  the  bubble  pulse  is  in  the  vicinity  of  one  of  the  natural 
periods  of  hull  flexure.  Indeed,  if  the  overall  dissipative  forces  (damping  forces)  are  low, 
whipping  response  magnitudes  (and  therefore  resulting  hull  damage)  could  be  quite 
substantial. 

To  see  potential  dynamic  effects  associated  with  matched  bubble  pulse  (i.e.  bubble 
pulse  frequency  matched  to  ship's  flexural  frequency),  explosive  charge  characteristics 
(charge  weight,  depth,  etc.)  must  first  be  defined.   Using  the  formulation  for  bubble 
pulsation  and  migration  discussed  in  chapter  3,  and  the  known  "wet"  natural  frequencies 
of  hull  flexure  (Table  4-2),  explosive  charge  characteristics  could  be  iteratively 
determined.  Table  4-3  specifies  explosive  charge  characteristics  used  in  this  analysis  to 
match  bubble  pulse  frequency  with  1st  and  2nd  horizontal  flexural  mode  natural 
frequencies  (i.e.  modes  6  and  8  in  Table  4-2).  A  charge  depth  of  100  feet  was  selected,  as 
deeper  depths  require  unrealistically  large  charges  to  match  the  1st  fundamental  flexural 
frequency  (1.37  Hz).   Ship  location  relative  to  the  charge  was  chosen  to  promote 
excitation  of  the  mode  (both  frequency  and  shape)  of  interest. 

Figure  4-4  shows  plots  of  bubble  radius,  depth,  and  migration  velocity  for  the 
explosive  charge  weight/depth  combination  selected  to  excite  mode  2  horizontal  whipping 
of  the  SSN.   It  can  be  seen  from  the  plot  of  radius  vs.  time  that  the  bubble  period  is  not 
constant  with  time,  but  increases  as  the  bubble  migrates  toward  the  surface    This  is  an 
important  phenomenon,  which  limits  the  ability  of  the  pulsating  bubble  to  excite  a 
particular  flexural  natural  frequency  for  a  prolonged  period  of  time  (thus,  the  bubble  pulse 
could  only  be  considered  to  be  quasi-periodic  in  nature). 

The  whipping  response  of  the  SSN  is  predicted  for  the  cases  when  excited  by 
explosion  bubble  pulse  characterized  by  charges  specified  in  Table  4-3.  Subroutine 
BUBBLE. M  predicts  explosion  bubble  radius  and  depth  time  histories  for  each  of  the 
charge  weight/depth  combinations  selected.   Subroutine  BUB  WHIP. M  predicts  the  nodal 
point  response  of  the  SSN,  initially  at  rest  and  subject  to  the  incident  bubble  pulse  loading. 
Appendix  D  presents  plots  of  displacements,  velocities,  and  accelerations  for  nodes  7,  10, 
and  20  for  the  charges  designed  to  excite  horizontal  modes  1  and  2  flexural  whipping. 
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Table  4-3.  Charge  characteristics  selected  to  excite  mode  1  and  mode  2  horizontal 


flexural  whipping  for  SSN. 

First  horizontal 
flexural  mode 

Second  horizontal 
flexural  mode 

Charge  weight  (lb  TNT) 

1000 

85 

Charge  depth  (ft) 

100 

100 

Bubble  pulse  frequency  (Hz) 

1.35 

3.05 

Depth  of  ship  (ft) 

70 

80 

Charge  location  along  length  of 
ship  (ft) 

157.5  (node  10) 

105.0  (node  7) 

Charge  standoff  from  ship 
centerline  (ft) 

33.5 

15 

Charge  depth  below  ship 
centerline  (ft) 

30 

20 

Because  of  the  swiftness  with  which  the  (large)  bubble  migrates  to  the  surface 
from  its  initial  depth  for  the  larger/lower  frequency  charge  (mode  1  whipping),  the 
calculated  response  time  period  was  limited  to  the  first  2.5  seconds.  Despite  the  short 
time  of  excitation  (only  3  flexure/bubble  periods),  it  can  be  easily  observed  that  the  quasi- 
periodic  bubble  pulse  does  indeed  magnify  the  whipping  response.  For  the  smaller/higher 
frequency  charge  (mode  2  whipping),  the  magnification  effect  can  be  more  easily  observed 
over  a  larger  number  of  whipping  cycles. 
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Figure  4-4.  Plots  of  bubble  radius,  depth,  and  migration  velocity  for  the  explosive  charge 
weight/depth  combination  selected  to  excite  mode  2  horizontal  whipping  of  the  SSN. 


90 


4.5        Damping  from  discrete  internal  structures  and  sloshing  liquids. 

Discrete  internal  structures.  As  discussed  in  chapters  2  and  3,  whipping  energy 
can  be  transferred  to  internal  structures  within  the  hull,  particularly  when  the  internal 
structure  is  "tuned"  to  the  bubble  pulse  frequency  (and  one  of  the  natural  frequencies  of 
hull  flexure).  Using  the  model  developed  in  chapter  3  (which  models  the  internal  structure 
as  a  discrete  mass  or  dynamic  vibration  absorber),  a  brief  analysis  is  conducted  to  gain 
some  insight  into  potential  dynamic  effects  of  these  internal  structures  on  the  whipping 
response  of  a  hull. 

The  analysis  discussed  in  the  previous  section  is  extended  to  include  "tuned" 
vibration  absorbers  (internal  masses  on  springs)  where  each  absorber  is  "tuned"  to  the 
matched  bubble  pulse/whipping  frequency    The  effects  of  a  "tuned"  internal  absorber  are 
predicted  for  the  case  of  the  bubble  pulse  matched  to  the  2nd  horizontal  flexural 
frequency.  The  internal  absorber  is  located  at  node  7  (i.e.  "attached"  to  node  7)  and  free 
to  move  in  the  horizontal  direction  only,  as  discussed  in  section  3.3.  Thus,  the  effect  of 
the  absorber  is  on  the  2nd  horizontal  whipping  mode 

The  effect  of  a  200  long  ton  internal  mass,  "tuned"  to  the  2nd  horizontal  whipping 
mode  and  "attached"  to  the  hull  at  node  7,  is  considered.  Appendix  E  presents  plots  of 
horizontal  displacements,  velocities,  and  accelerations  for  nodes  7,  10,  and  20,  as  well  as 
horizontal  displacement,  velocity,  and  acceleration  of  the  absorber  mass.  Figure  4-5  plots 
horizontal  displacement  of  node  7,  with  and  without  the  200  ton  internal  mass. 

Several  interesting  phenomenon  are  evident  from  the  plot  and  should  be  briefly 
noted.  First,  the  presence  of  the  absorber  has  a  significant  overall  effect  on  reducing  the 
magnitude  of  response.  While  no  specific  measure  of  "damping"  can  be  made,  there  is 
obviously  significant  energy  transferred  into  the  absorber.   Second,  the  ability  of  the 
absorber  to  "absorb"  whipping  energy  is  very  sensitive  to  the  frequency  of  excitation.  As 
discussed  in  the  previous  section,  the  bubble  period  (and  therefore  the  frequency  of 
whipping  excitation)  is  not  constant,  but  rather  increases  as  the  bubble  migrates  toward 
the  surface.  By  "tuning"  the  internal  absorber  to  the  2nd  horizontal  flexural  frequency,  the 
absorber  is  not  then  truly  tuned  to  the  excitation  frequency  (recall  that  matching  bubble 
period  to  hull  flexural  frequency  was  accomplished  by  iteration  of  average  bubble  periods 
-  over  the  3  seconds  of  bubble  pulsation).  Thus,  the  usual  "internal  absorber  effect"  (i.e. 
reducing  vibration  amplitude  of  the  main  structure  to  zero)  would  not  occur.  Figure  4-6 
plots  the  horizontal  displacement  of  node  7  and  the  horizontal  displacement  of  the 
absorber  mass.  Because  the  oscillation  of  node  7  and  the  oscillation  of  the  absorber  mass 
occur  at  slightly  different  frequencies,  a  sharp  (but  temporary)  reduction  in  node  7 
amplitude  can  be  seen  around  2.5  seconds,  when  the  absorber  shifts  from  a  lag  angle  to  a 
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lead  angle  (i.e  the  node  and  the  mass  are  temporarily  180°  out  of  phase)    What  can  be 
deduced  from  this  is  that  internal  structures  do  have  potential  for  significant  absorption  of 
hull  whipping  energy,  even  though  bubble  pulse  frequencies  are  not  truly  constant. 
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Figure  4-5    Horizontal  displacement  of  node  7,  with  and  without  internal  absorber 
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Figure  4-6.  Horizontal  displacement  of  node  7  and  200  LT  absorber  mass. 
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Liquid  sloshing.  As  discussed  in  chapters  2  and  3,  the  fundamental  mechanism  by 
which  energy  is  "transferred"  to  internal  liquids  through  the  mechanism  of  sloshing  can  be 
thought  of  much  in  the  same  way  as  that  by  which  it  is  transferred  to  internal  structures 
In  fact,  liquid  sloshing  may  be  modeled  quite  accurately  using  a  mechanical  equivalent 
system  approach,  where  each  sloshing  mode  is  considered  to  have  its  own  modal  mass, 
damping  and  stiffness,  each  responding  as  its  own  dynamic  system  (i.e.  its  own  dvnamic 
absorber).  Using  the  formulation  for  natural  sloshing  frequency,  modal  mass,  and  modal 
stiffness  presented  in  chapter  3  for  lateral  sloshing  in  a  rectangular  tank,  the  discussion  of 
damping  through  liquid  sloshing  follows  quite  nicely  behind  that  of  the  discrete  internal 
structure.  Recalling  equations  (3-35)  through  (3-37),  it  is  clear  that  the  natural  frequency, 
modal  mass,  and  modal  stiffness  of  each  mode  of  sloshing  are  dependent  only  upon  the 
mode  number  (n),  the  total  mass  of  liquid  in  the  tank  (M,),  and  the  geometry  of  the  tank 
(i.e.  tank  depth  and  width) 

Using  the  standard  rectangular  tank  formulations,  it  is  a  trivial  procedure  to 
generate  plots  of  natural  frequency  vs.  sloshing  mode  and  sloshing  modal  mass  vs. 
sloshing  mode.  When  this  is  done  over  a  wide  range  of  possible  tank  configurations,  it 
quickly  becomes  clear  that  liquid  sloshing  could  not  become  a  sigtnficant  contributor  to 
the  dissipation  or  absorption  of  hull  whipping  energy.  Figure  4-7  plots  both  natural 
sloshing  frequency  (con)  and  sloshing  modal  mass  (mn)  for  a  20'  x  20'  x  20'  rectangular 
tank  (a  total  liquid  weight  of  228.5  long  tons).  Two  interesting  things  may  be  noted  for 
this  very  large  tank  (much  larger  than  any  single,  unbaffled  tank  on  any  submarine).  First, 
the  only  sloshing  modal  mass  of  any  substance  (i.e.  which  might  effect  hull  whipping)  is 
that  of  the  first/fundamental  sloshing  mode  (i.e.  60  long  ton).   Second,  the  natural 
frequency  of  this  fundamental  sloshing  mode  is  only  2.24  rad/sec  (0.35  Hz),  well  below 
any  natural  flexural  frequency    As  other  tanks  are  considered,  it  is  found  that  the 
fundamental  sloshing  modes  (of  all  reasonable  tank  sizes)  fall  well  below  1  Hz.  While 
higher  sloshing  modes  may  have  natural  sloshing  frequencies  in  the  range  of  hull  whipping 
frequencies,  the  modal  masses  of  these  modes  are  significantly  smaller  than  would 
contribute  to  any  noticeable  reduction  in  whipping  response.   Appendix  F  presents  plots  of 
natural  sloshing  frequency  vs.  sloshing  mode  and  sloshing  modal  mass  vs.  sloshing  mode 
for  several  representative  rectangular  tank  geometries. 

From  this  basic  analysis  of  natural  sloshing  frequencies  and  sloshing  modal  masses, 
it  can  be  justifiably  deduced  that  liquid  sloshing  would  have  little  impact  on  the  whipping 
response  of  most  ships.  Because  the  significant  sloshing  modal  masses  have  natural 
frequencies  well  below  hull  flexural  frequencies,  it  can  also  be  deduced  that  the  major 


effect  of  liquids  with  free  surface  would  only  be  toward  the  addition  of  total  ship  mass  (i.e. 
the  significant  sloshing  modes  would  act  only  as  rigid  bodies) 
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4.6        Comparison  of  computational  algorithm  with  the  "Red  Snapper"  model  test. 

The  computational  algorithm  for  submarine  whipping  developed  in  chapter  3,  and 
applied  to  the  whipping  of  a  notional  SSN  in  the  previous  sections,  is  used  to  predict  the 
response  of  the  U.S.  Navy  test  platform  "Red  Snapper".  Thus,  some  "verification"  of  the 
computational  algorithm  is  made  with  actual  model  test  data.  Due  to  the  security 
concerns  with  the  actual  test  data  used  for  comparison,  the  following  is  presented  in  non- 
specific response  units  (i.e.  the  y-axes  are  unlabeled) 
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The  "Red  Snapper"  model  is  on  the  order  of  1/4  the  size  of  the  notional  SSN 
discussed  previously  (geometrically),  and  thus  natural  frequencies  of  hull  flexure  are  much 
higher  than  the  notional  SSN.  Whipping  response  was  measured  in  terms  of  local  axial 
strains  (at  specific  beam  element  locations)  and  horizontal  and  vertical  velocities  (at 
specific  node  locations)    For  comparative  purposes,  subroutine  BUB  WHIP  M  includes  a 
scheme  for  calculating  local  axial  strains  at  specified  finite  beam  elements    This  scheme  is 
based  upon  simple  beam  theory  an  provides  local  axial  strains  at  port  and  starboard 
extreme  fibers,  as  well  as  crown  and  keel  extreme  fibers  for  the  specified  elements  (see 
subroutine  BUB  WHIP.  M  in  Appendix  A) 

Three  specific  model  tests,  which  were  conducted  on  the  "Red  Snapper"  model, 
are  used  for  comparison    One  test  was  designed  such  that  the  explosion  bubble  pulse 
would  excite  mode  1  horizontal  whipping.  A  second  test  was  designed  such  that  the 
bubble  pulse  would  excite  mode  1  whipping  in  both  horizontal  and  vertical  planes.  A  third 
test  was  designed  such  that  the  bubble  pulse  would  excite  mode  2  horizontal  whipping. 
The  design  of  the  bubble  pulse  to  match  the  flexural  frequencies  and  mode  shapes  was 
carried  out  in  a  manner  similar  to  that  discussed  in  section  4.4.   In  all  cases,  there  where 
no  internal  structures  or  liquids,  and  thus  all  forces  were  either  hydrodynamic  or  structural 
in  nature 

General  comparisons.  Figure  4-8  shows  starboard-side  fiber  strain  histories  for 
elements  6  and  9  for  Test  1  (horizontal  mode  1  whipping),  and  strains  calculated  using  the 
computational  algorithm.  In  general,  it  appears  that  calculated  strains  match  well  with 
measured  through  the  2nd  bubble  period,  but  then  over  predict  strains.  Figure  4-9  shows 
horizontal  velocity  histories  for  nodes  10  and  20  for  Test  1,  and  velocities  calculated  using 
the  computational  algorithm.  Calculated  velocities  for  this  test  are  only  slightly  under 
predicted. 

Comparisons  between  measured  and  calculated  strains  and  velocities  for  Test  2 
(mode  1  horizontal  and  vertical  whipping)  and  Test  3  (mode  2  horizontal  whipping)  are 
presented  in  Appendix  G.   In  general,  the  computational  algorithm  slightly  under  predicts 
and  predicts  slightly  early  for  both  strains  and  velocities  for  Test  2.  For  Test  3,  the 
computational  algorithm  predicts  quite  well  through  the  first  two  bubble  periods,  but  then 
grossly  over  predicts  both  strains  and  velocities. 
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Figure  4-8.   Starboard-side  fiber  strain  histories  for  elements  6  and  9  for  Test  1  of 
the  "Red  Snapper"  model  test  (horizontal  mode  1  whipping),  and  strains  calculated  using 

the  computational  algorithm 


96 


Horizontal  velocity  vs.  time,  node  10 
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Figure  4-9.  Horizontal  velocity  histories  for  nodes  10  and  20  for  Test  1  of 
the  "Red  Snapper"  model  test  (horizontal  mode  1  whipping),  and  velocities  calculated 

using  the  computational  algorithm 
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Discussion  of  general  comparison.  There  are  several  possible  explanations  for  the 
failure  of  the  computational  algorithm  to  properly  predict  strains  and  velocities  for  the 
"Red  Snapper"  tests 

a.  The  bubble  model  used  for  formulation  of  equations  governing  bubble  dynamics 
has  been  known  to  be  inaccurate  past  the  2nd  bubble  period97.  Although  no  satisfactory 
purely  physical  model  exists  for  long-time  bubble  dynamics  (some  more  restrictive 
empirical  relations  are  available),  this  particular  bubble  model  has  proven  to  be  fairly 
accurate  through  the  2nd  bubble  period    The  reasons  for  the  inaccuracy  past  the  2nd 
bubble  period  lie  in  the  inability  of  the  model  to  properly  account  for  dissipation  of  bubble 
energy  during  bubble  pulsation  and  migration.  Thus,  the  model  believes  that  the  bubble 
maintains  more  of  its  original  energy  throughout  the  longer  time  history    This  translates 
into  greater  bubble  radii  and  greater  velocities  being  modeled  in  the  computational 
algorithm  (resulting  in  over  predicting  hydrodynamic  loading  past  the  2nd  bubble  period). 
Thus,  some  explanation  for  overproduction  of  response  past  the  2nd  bubble  period  exists. 

b    The  use  of  the  hydrodynamic  "line-structure"  model  in  the  computational 
algorithm  (i.e.  hydrodynamic  loading  being  applied  at  the  nodes  along  the  ship  centerline, 
vice  on  the  outer  hull  of  the  ship),  leads  to  inaccuracies  in  hydrodynamic  loading.  This  is 
particularly  important  when  the  explosive  charge  is  close  aboard    In  this  case,  the  actual 
hydrodynamic  loading  is  more  appropriately  seen  directly  on  the  outer  hull,  while  the 
computational  algorithm  believes  the  loading  to  be  applied  at  the  ship's  centerline  (i.e. 
calculated  "free-field"  fluid  velocities  and  accelerations).  Thus,  an  error  is  introduced  due 
to  the  unaccounted-for  distance  between  the  centerline  and  the  outer  hull  of  the  ship    For 
large  ships,  subjected  to  charges  with  large  standoffs,  this  geometric  inaccuracy  is  less 
important.  But,  for  smaller  ships  subjected  to  near-standoff  charges,  this  inaccuracy  could 
lead  to  significant  error. 

c.  There  is  a  general  inability  of  Morison's  equation  to  adequately  account  for  the 
transient-type  vortex  shedding  and  convection  (viscous  drag  forces),  which  is  undoubtedly 
important  in  this  application.  The  original  Morison's  equation  was  designed  to  account  for 
hydrodynamic  loading  for  the  relatively  low-frequency  wave  loading  on  pile  structures  in 
regular  sea  waves  (i.e.  quasi-steady)    There  is  no  adequate  data  in  the  literature  to 
support  extension  of  the  Morison  formulation  to  transient-type  analyses.  This  is 
particularly  important  in  the  justification  and  selection  of  the  coefficients  used  in  Morison's 
equation.  As  discussed  in  chapter  3,  coefficients  for  the  computational  algorithm  were 
merely  extrapolated  from  low  frequency  data.  Additionally,  this  data  was  itself  obtained 


97Misovec.  A..  "Comments  on  SITE  meeting  of  5/10/94".  (Unpublished  meeting  minutes).  1994.  pp.  1-2. 
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by  Fourier  time-averaging  of  measured  forces  (i.e.  quasi-steady).  For  transient-type 
analyses,  such  as  submarine  whipping,  a  more  complex  approach  would  probably  be 
appropriate  to  properly  account  for  the  transient  vortex  shedding  and  convection. 
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5.  CONCLUSIONS 

An  analytical  approach  has  been  presented  which  predicts  the  elastic  whipping 
response  of  a  submerged  submarine  to  the  pulsating  bubble  of  a  nearby  underwater 
explosion.  The  approach  presented  concentrates  on  mechanisms  associated  with  energy 
dissipation  (damping),  and  includes  specific  mechanisms  associated  with  hull  damping 
(material  and  structural),  internal  damping  (internal  structures  and  sloshing  liquids),  and 
external  damping  (hydrodynamic  damping  forces  in  conjunction  with  traditional  inertial 
loading  forces).  The  formulation  utilizes  a  numerical  time-domain  solution  technique  to 
directly  solve  the  nonlinear  equations  of  motion. 

While  the  general  mechanisms  associated  with  damping  are  reasonably  well 
understood,  application  of  this  understanding  to  complex  high-frequency,  non-harmonic 
vibration  problems  (such  as  ship  whipping  due  to  explosion  bubble  pulse  loading)  has  not 
been  met  with  much  success.  For  this  reason,  the  approach  has  been  designed  to  utilize  a 
robust  solution  technique  (the  time-domain  solution),  and  incorporate  physically-based 
models  where  possible. 

Modal  damping  ratios  for  material  and  structural  damping  were  assumed  to  be  on 
the  order  of  1  %  of  critical  for  fundamental  flexural  modes,  based  upon  data  available 
from  steel  marine  structures.  Thus,  results  for  "dry"  whipping  only  reflect  this  initial 
assumption.  With  this  initial  assumption  for  material  and  structural  damping,  "equivalent" 
modal  damping  ratios  for  mode  1  and  mode  2  whipping  (the  first  two  fundamental  flexural 
modes)  for  the  case  where  hydrodynamic  damping  effects  are  included  (via  Morison's 
equation),  are  computed  to  be  on  the  order  of  2.2  %. 

The  effects  of  internal  structures  and  sloshing  liquids  on  whipping  response  is 
approximated  using  a  200  long  ton  "tuned"  mass  damper,  located  to  effect  mode  2 
whipping.  The  effects,  illustrated  in  figure  4-5,  are  significant  only  when  bubble  pulse 
frequencies  are  precisely  equal  to  the  "tuned"  frequency  of  the  internal  structure.    The 
effects  of  internal  sloshing  liquids  are  shown  to  be  negligible  because  of  the  low  natural 
fundamental  sloshing  frequency  associated  with  larger  liquid  tanks. 

Several  of  the  selected  modeling  techniques  are  shown  by  comparison  with  actual 
test  data  to  be  questionable.  Firstly,  the  representation  of  hydrodynamic  forces  solely  by  a 
Morison's  equation  approach  fails  to  fully  account  for  the  expected  vortex  shedding  and 
convection  forces  during  the  transient  vibration.  The  selection  procedure  (or  more 
appropriately,  the  calculation  procedure)  of  hydrodynamic  coefficients  appropriate  for  ship 
whipping  must  be  better  understood  and  implemented.  Mere  extrapolation  of  known  data 
based  upon  harmonic,  low  frequency  motion  of  small  scale  cylinders  and  flat  plates  is 
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probably  inaccurate.  A  more  robust  model  for  hydrodynamic  forces,  coupling  numerical 
simulation  with  experimental  verification  should  be  developed. 

A  second  modeling  technique  of  questionable  success  is  the  bubble  pulsation  and 
migration  equations  used  to  produce  the  free-field  fluid  velocities  and  accelerations. 
Because  the  model  does  not  account  for  proper  dissipation  of  energy  by  the  pulsating 
explosion  bubble,  it  is  believed  that  the  model  over  predicts  fluid  loading  after  the  2nd 
bubble  period.  A  more  robust  model  for  fluid  loading  from  the  pulsation  of  the  explosion 
bubble,  again,  coupling  numerical  simulation  with  experimental  verification  should  be 
more  successful. 

Finally,  the  general  nature  of  the  line-structure  model  for  the  ship  structural  and 
hydrodynamic  interaction  (i.e.  the  fluid-structure  interaction)  introduces  geometric  errors, 
particularly  when  the  explosion  is  close  aboard.  A  more  robust,  3-dimensional  fluid- 
structure  interaction  solution  would  be  more  precise  for  the  solution  process. 

Even  with  several  questionable  physical  models  for  particular  aspects  of  the 
whipping  model,  comparison  with  test  data  shows  that  the  model  provides  reasonable 
results,  particularly  in  early-time  response.  It  is  believed  that  the  general  approach  of 
solving  the  equations  of  motion  directly  in  the  time  domain  was  reasonable  and  necessary. 
It  is  in  the  modeling  of  each  individual  component  of  the  forces  on  the  whipping  ship 
where  future  work  should  be  concentrated. 
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%  Subroutine  IN.M 

% 

%  This  matlab  subroutine  is  written  as  an  INPUT  FILE  for  all  other  MATLAB 

%  subroutines.  Data  is  loaded  from  formatted  batch-type  input  files  using 

%  the  MATLAB  "load"  command,  and  arrays  and  constants  are  generated  for 

%  the  other  MATLAB  subroutines. 

% 

%  MISCELLANEOUS/DEFAULT  VALUES 

% 

%  Number  of  hull  stations  or  finite  beam  elements 

NE=19; 

%  Number  of  hull  sections  or  finite  beam  nodes 

NNODES=20: 

%  Number  of  degrees  of  freedom  (4  DOF  at  each  node) 

NDOF=4*NNODES; 

% 

% 

%  PROPERTIES  OF  EACH  FINITE  BEAM  ELEMENT 

load  element.dat 

% 

%  Element  number 

e=element(l.:): 

%  Cross-sectional  area  of  element  (in2) 

Ae=element(2.:); 

%  Length  of  element  (ft) 

Le=element(3.:); 

%  Area  effective  in  y-directed  shear  (in2) 

Asy=element(4.:); 

%  Area  effective  in  z-directed  shear  (in2) 

Asz=element(5.:); 

%  Second  moment  of  area  about  y-axis  (in2ft2) 

Iay=element(6.:); 

%  Second  moment  of  area  about  z-axis  (in2ft2) 

Iaz=element(7.:); 

%  Y-distance  from  the  neutral  axis  to  the  outer  fiber,  of  the  crowti  (ft) 

Yc=element(8.:): 

%  Y-distance  from  the  neutral  axis  to  the  outer  fiber  of  the  keel  (ft) 

Yk=element(9.:); 

%  Z-distance  from  the  centerline  to  the  outer  fiber  of  the  port  hull  (ft) 

Zp=element(10.:); 

%  Z-distance  from  the  centerline  to  the  outer  fiber  of  the  stbd  hull  (ft) 

Zs=element(ll.:); 

% 

% 

%  PROPERTIES  OF  EACH  HULL  SECTION  (NODE) 

load  node.dat 

% 

%  Node  number 

nnode=node(l.:): 

%  Length  of  each  hull  section  (ft) 

Lh=node(2.:); 

%  Weight  of  hull  section  (i.e.  in  air)  (lton) 

wh=node(3.:): 

%  mass  of  hull  section  (lbf-sec2/ft=slug) 
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mh=wh*2240./32.174; 

%  Buoyancy  of  hull  section  (displacement  of  salt  water)  (lton) 

bh=node(4.:): 

%  Weight  moment  of  inertia  of  hull  section  about  y-axis  (lton-ft2) 

Imyh=node(5.:): 

%  Mass  moment  of  inertia  of  hull  section  about  y-axis  (slug-ft2) 

Imy=Irnyh*2240./32.174; 

%  Weight  moment  of  inertia  of  hull  section  about  z-axis  (lton-ft2) 

Imzh=node(6.:); 

%  Mass  moment  of  inertia  of  hull  section  about  z-axis  (slug-ft2) 

Imz=Imzh*2240./32.174; 

%  Maximum  section  beam  of  main  hull  perp.  to  motion  in  z  direction  (ft) 

Bz=node(7.:); 

%  Maximum  section  beam  of  main  hull  perp.  to  motion  in  y  direction  (ft) 

By=node(8.:): 

%  Cross-sectional  area  of  Appended  hull  (appended  portion  of  each  section) 

%  perp.  to  motion  in  y  direction  (ft2) 

Aapy=node(9.:); 

%  Cross-sectional  area  of  Appended  hull  (appended  portion  of  each  section) 

%  perp.  to  motion  in  z  direction  (ft2) 

Aapz=node(10.:); 

%  Width  of  Appended  hull  perp.  to  motion  in  y  direction  (ft) 

Dap\-node(ll.:): 

%  Width  of  Appended  hull  perp.  to  motion  in  z  direction  (ft) 

Dapz=node(12.:); 

%  Length  of  Appended  hull  perp.  to  motion  in  y  direction  (ft) 

Lapy=node(13,:); 

%  Length  of  Appended  hull  perp.  to  motion  in  z  direction  (ft) 

Lapz=node(14.:): 

% 

%  Hvdrodvnamic  coefficients  (non-dimensional) 

% 

%  Added  mass,  inertia,  and  drag  coefficients  for  each  main  hull  section 

%  for  motion  in  vertical  (y)  and  horizontal  (z)  directions 

Caly=node(15.:); 

Calz=node(16.:); 

Cmly=node(17.:); 

Cmlz=node(18.:); 

Cdly=node(19.:): 

Cdlz=node(20.:); 

%  Added  mass,  inertia,  and  drag  coefficients  for  appended  hull  sections 

%  (appended  portion  of  each  hull  section)  for  motion  in  y  and  z  directions 

Ca2y=node(21.:); 

Ca2z=node(22.:): 

Cm2y=node(23.:); 

Cm2z=node(24,:); 

Cd2y=node(25.:); 

Cd2z=node(26.:); 

% 

% 

%  MISCELLANEOUS  INPUT 

load  misc.dat 

% 

%  Acceleration  due  to  gravity  (ft/sec2) 
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g=misc(l): 

%  Density  of  salt  water  (Ibf-sec2/ft4=slugs/ft3) 

rho=misc(2): 

%  Young's  modulus  for  steel  (lbf/in2) 

E=misc(3); 

%  Poisson  ratio  for  steel  (non-dimensional) 

Nu=misc(4): 

%  Yield  stress  for  steel  (lbf/in2) 

Sigy=misc(5): 

%  Shear  modulus  for  steel  (lbf/in2) 

G=misc(6); 

%  Assumed  material  (hysteresis)  modal  viscous  damping  coefficients  for 

%  all  low  frequency  (fundamental)  flexural  "dry"  modes  (non-dimensional) 

damp=misc(7); 

%  Whipping  mode  to  scale  as  initial  condition  for  use  in  FREEWHTP.M 

MODE=misc(8); 

%  Switch  to  determine  whether  or  not  "wet"  whipping  is  used  in  FREEWHTP.M 

WET=misc(9); 

%  Node  numbers  to  plot  displacements.velocities.accelerations 

NODEl=misc(10); 

NODE2=misc(ll); 

NODE3=misc(12); 

NODE4=misc(13); 

%  Element  numbers  to  plot  strains 

ELEMENT  l=misc(  14); 

ELEMENT2=misc(15); 

ELEMENT3=misc(16); 

% 

%  Input  parameters  for  Newmark  time  integration  algorithms 

% 

%  time  step  (sec) 

dt=misc(17); 

%  gamma  parameter  (non-dimensional) 

gammap=misc(18); 

%  beta  parameter  (non-dimensional) 

betap=nusc(19); 

%  time  limit  (sec) 

tma.\=misc(20); 

%  vector  norm  tolerance 

tol=misc(21): 

% 

%  Input  for  internal  vibration  absorbers 

% 

%  "Reasonable"  weight  of  an  absorber  or  internal  structure  (lton) 

wa=misc(22): 

%  mass  of  absorber  or  internal  structure  (slug) 

ma=wa*2240/32.174; 

%  Node  number  where  absorber  is  "attached"  to  the  hull 

nabs=misc(23): 

% 

%  Input  for  internal  sloshing  tank 

% 

%  "Reasonable"  weight  of  liquid  in  a  (large)  internal  tank  (lton) 

wl=rrusc(24): 
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%  mass  of  liquid  (slug) 

Ml=wl*2240./32.174; 

%  "Reasonable"  aspect  ratio  of  tank  depth/width  (h/a)  (non-dimensional) 

ha=misc(25): 

%  Node  number  where  tank  is  "attached"  to  the  hull 

ntank=misc(26); 

% 

%  Input  for  explosion  bubble  dynamics 

% 

%  Initial  depth  of  charge  (ft) 

Dchg=misc(27); 

%  Charge  weight  of  TNT  (lb) 

Wchg=misc(28): 

%  Initial  location  of  charge  relative  to  ship  hull  (in  hull  coordinate  system) 

xchg=misc(29); 

ychg=misc(30); 

zchg=misc(31); 

%  Initial  depth  of  the  ship  (ft) 

Dship=misc(32); 
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%  SUBROUTINE  HULL.M 

% 

%  This  subroutine  performs  the  following: 

%  1.  Calculates  the  global  mass  and  stiffness  matrices  (K  and  M) 

%  for  the  3-D  beam  finite  element  formulation  (ignoring  axial 

%  and  torsional  effects) 

%  2.  Solves  the  generalized  eigenvalue  problem  to  obtain  the 

%  "dr>  mode"  eigenvalues  (natural  freqs  squared)  and 

%  corresponding  eigenvectors  ("dry  mode"  shapes). 

%  3.  Calculates  the  Rayleigh  coefficients  for  critical  hull 

%  (structural  and  hull  hysteresis)  damping  defined  by 

%  the  1st  and  3rd  lowest  fundamental  flexural  modes 

%  4.  Calculates  the  equivalent  viscous  hull  clamping  matrix  for  the 

%  calculated  Rayleigh  coefficients. 

%  5.  Calculates  vectors  for  modal  mass,  modal  damping,  and 

%  modal  stiffness  (for  N  modes)  corresponding  to  the  "dry 

%  mode"  shape  vectors    These  are  used  for  optimization 

%  of  internal  absorber  and  liquid  sloshing  dampers  (tanks) 

%  in  other  subroutines. 

%  6.  Provides  for  validation  against  the  "exact"  solution  given  by 

%  Rao  (chapter  8)  for  a  continuous  uniform  free-free  beam. 

%  The  natural  frequencies  of  the  finite  element  solution 

%  can  be  compared  to  the  natural  frequencies  of  the 

%  "exact"  solution.  This  may  be  disabled  after  validation 

%  bv  placing  %  before  each  calculation. 

% 

%(D 

%  Formulate  the  global  (hull)  stiffness  and  mass  matrices  (K  and  M). 

% 

%  Initialize  K  and  M 

K=zeros(NDOF); 

M=zeros(NDOF); 

% 

%  Determine  the  global  stiffness  matrix  (Ibf/ft  for  translation  or  or  lbf-ft  for  rotation) 

% 

%  Determine  the  element  stiffness  matrices  (KE) 

% 

%  Step  through  each  of  the  19  elastic  beam  elements 

fore=l:NE 

% 

FV=(12*E*Iaz(e))/(G*Asy(e)*Le(e)A2); 

Pz=(12*E*Iay(e))/(G*Asz(e)*Le(e)A2); 

bz=E*Iaz(e)/((  l+Py)*Le(e)A3); 

by=E*Iay(e)/((  l+Pz)*Le(e)A3): 

Kl  l=[12*bz  0  0  6*bz*Le(e);  0  12*by  -6*by*Le(e)  0; 

0  -6*by*Le(e)  (4+Pz)*by*Le(e)A2  0;  6*bz*Le(e)  0  0  (4+Py)*bz*Le(e)A2]; 

K21=[-12*bz  0  0  -6*bz*Le(e).  0  -12*by  6*by*Le(e)  0; 

0  -6*by  *Le(e)  (2-Pz)*b\ *Le(e)A2  0;  6*bz*Le(e)  0  0  (2-Py)*bz*Le(e)A2]; 

K12=K21\ 

K22=[12*bz  0  0  -6*bz*Le(e):  0  12*by  6*by*Le(e)  0; 

0  6*by*Le(e)  (4+Pz)*by*Le(e)A2  0;  -6*bz*Le(e)  0  0  (4+Py)*bz*Le(e)A2]; 

KE=[K11  K12;K21  K22]; 
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%  Transform  element  stiffness  matrices  into  their 

%  global  coordinate  matrices 

forie=l:8 

forje=l:8 

i=ie+4*(e-l); 
j=je+4*(e-l); 
K(ij)=K(ij)+KE(ieje); 
end 
end 
end 
% 

%  Once  each  of  the  elastic  beam  elements  has  been  considered  and  added. 
%  the  global  (hull)  stiffness  matrix  is  K. 
% 

%  Determine  the  global  mass  matrix  (Ibf-sec2/ft=slug  for  translation  or  slug-ft2  for  rotation) 
% 

%  Step  through  each  of  the  20  nodes 
% 
forn=l:NNODES 

%  Step  through  each  of  the  4  dof  at  each  node  and  define  the  NODE  mass  matrix 
MN(l.l)=mh(n); 
MN(2.2)=mh(n); 
MN(.3.3)=Imy(n); 
MN(4.4)=Imz(n); 
%  Transform  node  mass  matrices  for  each  node  into  their  global  coordinate  matrices 
for  in=l:4 

forjn=l:4 

i=in+4*(n-l); 
j=jn+4*(n-l); 
M(ij)=M(i.j)+MN(in.jn); 
end 
end 
end 
% 

%  Once  each  of  the  nodes  has  been  considered,  the  global 
%  mass  matrix  is  M 
% 

%(2) 

%  Solve  the  eigenvalue  problem  to  obtain  the  "dry  mode"  natural  frequencies  (square  roots  of 
%  eigenvalues)  and  corresponding  "drv  mode"  shape  vectors  (eigenvectors). 
% 

%  ***  It  should  be  noted  that  because  the  hull  has  4  rigid  body  modes. 
%     the  lowest  4  eigenvalues  will  be  zero  (zero  natural  frequency). 
%     These  correspond  to  the  lowest  4  eigenvectors  which  are  the 
%      RIGID  BODY  mode  shapes.  Thus,  the  lowest  FLEXURAL  mode 
%     would  be  mode  n=5. 
% 

%  "Phi  1"  is  a  (NxN)  matrix  w  hose  columns  are  "dry  mode"  shape  vectors. 

%  "Wl"  is  a  diagonal  (NxN)  matrix  whose  diagonal  elements  are  (square  of  the  natural  frequencies). 
%  (Note:  W  is  NOT  necessarilv  defined  in  ascending  order). 
% 

[Phil.Wl]=eig(K.M); 
% 
%  Define  a  VECTOR  of  natural  frequencies: 
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forr=l:NDOF 

wl(r)=sqrt(Wl(r.r)); 
end 
% 

%  Sort  natural  frequencies  and  define  a  sector  of  natural  frequencies  given  in  ascending  order: 
(w2.i]=sort(wl); 

%  where  "i"  is  a  vector  of  indices  such  that  w2=wl(i). 
% 

%  Sort  the  matrix  of  eigenvectors  (mode  shapes)  as  columns  in  an  order  corresponding  to  ascending 
%  natural  frequencies  (defined  by  index  vector  "i"): 
forr-l:NDOF 

s=i(r); 

Phi2(:,r)=Phil(:.s); 
end 
% 

%  Plot  the  2  rigid  body  mode  shapes  and  the  first  4  flexural  mode  shapes  corresponding  to  motion  in  the 
%  horizontal  plane  (for  example): 
%  Step  through  the  first  12  modes: 
%form=l:12 


% 

phi=Phi2(:.m); 

% 

forn=l:NNODES 

% 

nd(n)=n; 

% 

ny=l+4*(n-l): 

% 

nz=2+4*(n-l); 

% 

phiy(n)=phi(ny): 

% 

phiz(n)=phi(nz); 

% 

end 

% 

%  For  rigid  body  modes  (m<=4): 

% 

ifm<=4 

% 

figured ) 

% 

ifm<=2.  rbm=l; 

% 

elseifm<=4.  rbm=2; 

% 

end 

% 

ifphiz(l)~=0 

% 

subplot(2.1.rbm).  plot(nd.phiz) 

% 

axis([l  NNODES-1.0  1.0)) 

% 

title('Dry  mode  shape  of  rigid  body  mode  in  horizontal  (y=0)  plane' 

% 

end 

% 

%  For  flexural  modes  (m>4): 

% 

else 

% 

figure(2) 

% 

mcount=(m-4); 

% 

if  mcount<=2 

% 

flexmode=l; 

% 

elseif  mcount<=4 

% 

flexmode=2; 

% 

elseif  mcount<=6 

% 

flexmode=3; 

% 

else.    flexmode=4: 

% 

end 

% 

ifphiz(l)^0 

% 

subplot(4. l.flexmode).  plot(nd.phiz) 

% 

axis([l  NNODES-1.0  1.0]) 

% 

title(['Dry  mode  shape  of  flexural  mode  '.int2str(flexmode).'  in  ... 
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..horizontal  (y=0)  plane']) 
end 

end 
end 
Print  "dry"  mode  shapes 

figure(  1 ) 

set(gcf.'PaperPosition'.[1.0.  4.5.  6.5.  4.0]) 

print 

figure(2) 

set(gcf'PaperPosition'.[l  0.  2.0.  6.5.  8.5]) 

print 


(3) 


Calculate  the  Rayleigh  coefficients  defined  by  the  natural  frequencies  of  the  1st  and  3rd  flexural  modes 
(1st  and  2nd  flexural  mode  for  one  of  the  planes  of  motion  -  horizontal  or  vertical).  Ignoring  the  4 
rigid  body  modes,  these  would  be  correspond  to  modes  n=5  and  n=7. 


% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

%  Define  the  "dry  mode"  natural  frequencies  for  modes  5  and  7: 

w5=w2(5); 

w7=w2(7); 

% 

%  Find  the  Rayleigh  coefficients  using  equation  (3-16): 

w57=[w7  -\\5;  -l/w7  l/w5]; 

damping=[damp;  damp]; 

pre=2*(w5*w7)/(w7A2-w5A2); 

coefr=pre*Yv57*damping; 

alpha=coeff(l); 

beta=coeff(2); 

% 

% 

%  (4) 

%  Calculate  the  equivalent  viscous  damping  matrix  for  hull  damping  using  the  calculated  Rayleigh 

%  coefficients: 

% 

C=aJpha*M+beta*K; 

% 

% 

%(5) 

%  Define  the  "dry"  mode  shape  vectors  and  natural  frequencies  and  calculate  vectors  for  modal  mass 

%  (Mn),  modal  damping  (Cn).  and  modal  stiffness  (Kn)  corresponding  to  the  "dry  mode"  mode 

%  shape  vectors: 

% 

forn=l:NDOF 

phin=Phi2(:,n); 

phint=phin'; 

Phi(:.n)=phin; 

Mn(n)=phint*M*phin: 

Cn(n)=phint*C*phin; 

Kn(n)=phint*K*phin; 

w(n)=w2(n); 
end 
% 
%  (6) 


%  "dry"  mode  shape  vectors 
%  modal  mass  vector 
%  modal  damping  vector 
%  modal  stiffness  vector 
%  "dry"  natural  frequencies 
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%  Validation  of  a  sample  finite  element  formulation  (uniform  beam)  against  the  "exact"  solution,  for  a 

%  defined  axisvmmetnc  (about  x-axis)  beam  arrangement. 

% 

%  The  "exact"  solution  is  given  by  Rao  (Chapter  8)  for  a  uniform,  continuous 

%  free-free  beam  with  axisymmetrv  about  the  x-axis.  This  "exact"  solution 

%  does  not  include  shear  effects,  so  the  approximate  solution  should 

%  provide  natural  frequencies  which  are  slightly  lower  than  the  "exact" 

% 

%  Natural  frequencies  for  the  "exact"  solution  for  the  first  3  transverse  flexural  modes  (modes  1.2.3): 

% 

%  Bnl  1=4.73; 

%  Bnl2=7.85; 

%  Bnl3=11.0; 

%  rho=mh(l)/Lh(l): 

%  Sqroot=sqrt((E*Iay(l))/(rho*(Lh(l)*NE)A4)); 

%  omega  l=(BnllA2)*Sqroot; 

%  omega2=(Bnl2A2)*Sqroot: 

%  omega3=(Bnl3A2)*Sqroot; 

% 

%  Compare  "exact"  solution  to  the  (approximate)  numerical  solution. 

%  For  approximate  solution,  rigid  body  modes  are  ignored  and  the  1st.  2nd 

%  and  3rd  flexural  modes  for  transverse  displacement  (in  each  plane  (y=0,z=0)) 

%  are  considered.  It  should  be  noted  that  because  of  axisymmetrv.  natural 

%  frequencies  and  mode  shapes  in  each  plane  should  be  the  same 

%  (i.e.  repeated  eigenvalues). 

wl=w(5); 
w2=w(7); 
w3=w(9); 


% 

% 

% 

% 

%  Compare 

% 

omega  1 

% 

wl 

% 

omega2 

% 

v\2 

% 

omega3 

% 

w3 
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%  SUBROUTINE  WETMODE.M  (called  as  a  MATLAB  function) 

% 

%  This  subroutine  performs  the  following: 

%  1 .  Calculates  the  fluid  added  mass  matrix  for  the  submarine 

%  using  a  Monson  equation  approach  (i.e.  stnp  theory)  This 

%  approach  includes  added  mass  effects  in  the  translational 

%  degrees  of  freedom  only  (rotational  fluid  added  mass  effects 

%  are  not  included). 

%  2.  Adds  the  fluid  added  mass  matrix  to  the  hull  mass  matrix  (subroutine 

%  HULL.M)  and  solves  the  generalized  eigenvalue  problem  to  obtain 

%  the  "wet  mode"  eigenvalues  (natural  frequencies  squared)  and 

%  corresponding  "wet  mode"  eigenvectors  ("wet  mode"  shapes) 

%  3.  Calculates  vectors  for  "wet"  modal  mass,  modal  damping,  and  modal 

%  stiffness  (for  N  modes)  corresponding  to  the  "wet  mode"  shape 

%  vectors.  These  are  used  for  optimization  of  internal  absorber 

%  and  liquid  sloshing  dampers  (tanks)  in  other  subroutines 

% 

% 

%(D 

%  Calculate  the  fluid  added  mass  matrix  ( lbf-sec2/ft=slug  for  translation). 

% 

%  Initialize  the  fluid  added  mass  matrix 

Ma=zeros(NDOF); 

% 

%  Step  through  each  of  the  20  hull  sections  (nodes)  and  define  the  NODE  fluid  added  mass  matrices: 

forn=l:NNODES 

%  Calculate  the  fluid  added  mass  associated  with  the  main  hull  (Mai ) 
%  and  the  appended  hull  (Ma2) 
% 

%  Step  through  each  of  the  4  dof  at  each  node  and  define  the  NODE 
%  fluid  added  mass  matrices  for  the  main  hull  and  appended  hull 
MNal(l.l)=rho*(pi*By(n)A2/4)*Caly(n)*Lh(n); 
MNa2(l.l)=rho*Aapy(n)*Ca2y(n)*Lapy(n); 
MNal(2.2)=rho*(pi*Bz(n)A2/4)*Calz(n)*Lh(n). 
MNa2(2.2)=rho*Aapz(n)*Ca2z(n)*Lapz(n); 
MNal(3.3)=0; 
MNa2(3.3)=0; 
MNal(4.4)=0; 
MNa2(4.4)=0; 
%  Transform  NODE  fluid  added  matrices  into  their  global  coordinate  matrices 
for  in=l:4 

forjn=l:4 

i=in+4*(n-l); 
j=jn+4*(n-l); 

Ma(i.j)=Ma(i.j)+MNal(in.jn)+MNa2(injn); 
end 
end 
end 
% 

%  Once  each  of  the  nodes  has  been  considered,  the  global  fluid  added  mass  matrix  is  Ma 
% 

%  (2) 

%  Add  the  fluid  added  mass  matrix  to  the  hull  mass  matrix  and  solve  the  generalized  eigenvalue 
%  problem  to  obtain  the  "wet  mode"  eigenvalues  (natural  frequencies  squared)  and  corresponding 
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%  "wet  mode"  eigenvectors  ("wet  mode"  shapes). 

% 

%  Total  mass  matrix  (Mt): 

Mt=M+Ma; 

% 

%  "Phil"  is  a  (NxN)  matrix  whose  columns  are  "wet  mode"  shape  vectors. 

%  "Wl"  is  a  diagonal  (NxN)  matrix  whose  diagonal  elements  are  (square  of  the  "wet"  natural 

%  frequencies)  (Note:  W  is  NOT  necessarily  defined  in  ascending  order). 

% 

[Phil,Wl]=eig(K,Mt); 

% 

%  Define  a  VECTOR  of  natural  frequencies: 

forr=l:NDOF 

wl(r)=sqrt(WT(r.r)); 
end 
% 

%  Sort  natural  frequencies  and  define  a  vector  of  "wet"  natural  frequencies  given  in  ascending  order: 
[w2,i]=sort(wl); 

%  where  "i"  is  a  vector  of  indices  such  that  w2=wl(i). 
% 

%  Sort  the  matrix  of  eigenvectors  ("wet"  mode  shapes)  as  columns  in  an  order  corresponding  to  ascending 
%  "wet"  natural  frequencies  (defined  by  index  vector  "i"): 
forr=l:NDOF 

s=i(r); 

Phi2(:.r)=Phil(:.s); 
end 
% 

%  Plot  the  2  rigid  body  "wet"  mode  shapes  and  the  first  4  flexural  "wet"  mode  shapes  corresponding  to 
%  motion  in  the  horizontal  plane  (for  example): 
%  Step  through  the  first  12  modes: 
%form=l:12 


% 

phi=Phi2(:.m); 

% 

forn=l:NNODES 

% 

nd(n)=n; 

% 

ny=l+4*(n-l); 

% 

nz=2+4*(n-l): 

% 

phiy(n)=phi(ny); 

% 

phiz(n)=phi(nz); 

% 

end 

% 

%  For  rigid  body  modes  (m<=4): 

% 

% 

ifm<=4 

% 

figure(l) 

% 

if  m<=2,  rbm=l; 

% 

elseif  m<=4,  rbm=2; 

% 

end 

% 

ifany(phiz)~^0 

% 

subplot(2,l,rbm).  plot(nd.phiz) 

% 

axis([l  NNODES-1.0  1.0]) 

% 

title(' Wet  mode  shape  of  rigid  body  mode  in  horizontal  (y=0)  plane') 

% 

end 

% 

%  For  flexural  modes  (m>4): 

% 

else 

% 

figure(2) 
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% 

mcount=(m-4); 

% 

if  mcount<=2 

% 

fle.\mode=  1 : 

% 

elseif  mcount<=4 

% 

flexmode=2; 

% 

elseif  mcount<=6 

% 

flexmode=3; 

% 

else 

% 

flexmode=4; 

% 

end 

% 

ifphiz(l)~=<) 

% 

subplot(4.1.flexmode).  plot(nd.phiz) 

% 

axis([l  NNODES-I.O  l.OJ) 

% 

title(['Wet  mode  shape  of  flexural  mode  '.int2str(fle.\mode).'  in 

% 

...horizontal  (y==0)  plane' |) 

% 

end 

% 

end 

%  end 

%  Print  wet  mode  shapes 

% 

figure(  1 ) 

% 

set(gcf.'PaperPosition'.[1.0.  4.5.  6.5.  4.0]) 

% 

print 

% 

figure(2) 

% 

set(gcf.'PaperPosition'.[1.0.  2.0,  6.5.  8.5]) 

% 

print 

% 

%(3) 

%  Define  the  "wet"  mode  shape  vectors  and  natural  frequencies  the  first  and  calculate  vectors  for  "wet" 
%  modal  mass  (Mnw).  modal  damping  (Cnw).  and  modal  stiffness  (Knw)  corresponding  to  the  "wet" 
%  mode  shape  vectors  using  equations  (2-13): 


% 

for  n= 

=  l:NDOF 

phin=Phi2(:,n); 
phint^phin'; 

Phiw(:.n)=phin; 

% 

Mnw(n)=phint*Mt*phin; 

% 

Cnw(n)=phint*C*phin; 

% 

Knw(n)=phint*K*phin; 

% 

ww(n)=w2(n); 

% 

end 

% 

% 

% 

"wet"  mode  shape  vectors 
"wet"  modal  mass  vector 
"wet"  modal  damping  vector 
"wet"  modal  stiffness  vector 
"wet"  natural  frequencies 
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%  FREE  WHIP  M 

% 

%  This  subroutine  performs  the  following: 

%  Calculates  the  response  (nodal  point  displacements,  velocities. 

%  and  accelerations)  of  the  submarine  during  a 

%  free  vibration  with  initial  displacements  given  by  a 

%  scaled  mode  shape  vector  for  selected  mode  shapes 

%  Response  is  calculated  using  the  New  mark  time  integration 

%  method,  with  iteration  conducted  at  each  time  step  using  a 

%  Newion-Raphson  method 

%  Utilize  the  following: 

%  -  K.C.M.Phi  calculated  in  HULL.M 

%  -  Hydrodynamic  coefficients  provided  in  IN.M 

%  -  Parameters  for  Newmark  algorithm  provided  in  IN.M 

% 

%  Define  nodes  to  use  for  response  calculations 

N(l)=NODEl; 

N(2)=NODE2; 

N(3)=NODE3; 

N(4)=NODE4; 

%  Define  elements  to  use  for  axial  strain  calculations 

E11=ELEMENT1; 

E12=ELEMENT2; 

E13=ELEMENT3; 

% 

%  Set  initial  and  final  times  for  Newmark  integration,  and  limit  on  number 

%  of  Newton-Raphson  iterations  at  each  time  step 

t=0; 

tf=tmax. 

jlimit=50; 

% 

%  Initialize  the  loading  forces  on  the  hull 

%  Loading  forces  related  to  displacement 

Fdisp=zeros(l.NDOF)'; 
%  Loading  forces  related  to  velocity 

Fvel=zeros(l.NDOF)'; 
%  Loading  forces  related  to  acceleration 

Facc=zeros(  1  .NDOF)'; 
%  Calculate  (constant)  static  forces  (i.e.  hydrostatic/buoyancy  and  weight) 
Fgrav-zeros(l.NDOF)'; 
Fhs=zeros(l.NDOF)\ 
forn=l:NNODES 

ny=(n-l)*4+l; 

%  Hydrostatic/buoyancy  force  ( lbf) 

Fhs(ny)=bh(n)*2240; 

%  Static  force  of  gravity/weight  (lbf) 

Fgrav(ny)=-wh(n)*2240; 
end 

%  Total  loading  forces 
Floadt=Fdisp+Fvel+Facc+Fhs+Fgrav; 
% 

%  Define  initial  conditions  (displacements/velocities/accelerations) 
%  given  by  a  scaled  mode  shape  vector  (for  free  whipping) 
ifMODE-=l.  xt=(0.02)*sum(Lh)*Phi(:,6); 
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else.  xt=(0.02)*sum(Lh)*Phi(:.8); 

end 

% 

xdt=zeros(l.NDOF)'; 

Minv=inv(M); 

xddt=(Minv*(Floadt-C*\dt-K*xt)); 

% 

%  Perform  the  Newmark  time  integration 

% 

%  Step  through  time 

% 

%  Define  constants  to  be  used  in  Neuton-Raphson  iterations 

cl=(l/(betap*dtA2)h 

c2=(l/(betap*dt)); 

c3=(l/(2*betap)): 

c4=(gammap/(betap*dt)); 

c5=gammap/betap; 

c6=(  l-gammap/(2*betap)); 

c7=(l/(2*betap)-l); 

% 

d\=zeros(l.NDOF)'; 

% 

for  1=1:300; 

t=t+dt 

%  Set  initial  displacments.velocities.accelerations.  and  loading  forces  for 
%  the  current  time  step  (equal  to  the  final  for  the  last  time  step) 
xj=xt; 
xdj=xdt: 
xddj=xddt; 
Floadj=Floadt; 
% 

%  Perform  Newion-Raphson  iterations  at  current  time  step 
forj=l:jlimit 
% 

%  Calculate  loading  forces  for  the  current  iteration 
forn=l:NNODES 

%  Loading  forces  related  to  velocity 
ny=4*(n-l)+l; 
nz=4*(n-l)+2; 

Fvel(ny)=-0.5*rho*By(n)*Cdly(n)*xdj(ny)*abs(xdj(n\  ))*... 
...Lh(n)-0.5*rho*Dapy(n)*Cd2y(n)*xdj(ny)*abs(xdj(ny))*Lapy(n); 
Fvel(nz)=-0.5*rho*Bz(n)*Cdlz(n)*xdj(nz)*abs(xdj(nz))*... 
...Lh(n)-0.5*rho*Dapz(n)*Cd2z(n)*xdj(nz)*abs(xdj(nz))*Lapz(n); 
%  Loading  forces  related  to  acceleration 

Facc(ny)=-rho*(pi*By(n)A2/4)*Caly(n)*Lh(n)*... 
...xddj(ny)-rho*Aapy(n)*Ca2y(n)*Lapy(n)*xddj(ny); 
Facc(nz)=-rho*(pi*Bz(n)A2/4)*Calz(n)*Lh(n)*... 
...xddj(nz)-rho*Aapz(n)*Ca2z(n)*Lapz(n)*xddj(nz); 
end 
% 

%  Total  loading  force  for  the  current  iteration 
Floadj=Fdisp+F\el+Facc+Fhs+Fgrav; 
% 
%  Calculate  Jacobian  matrices  for  each  of  the  loading  forces 
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%  (square  diagonal  matrices  for  this  formulation) 

dFdxdj  1  =zeros(  1  .NDOF)': 

dFdxddj  l=zeros(  1  .NDOF)'; 

dFdxj=zeros(NDOF); 

dFd\dj=zeros(NDOF); 

dFd.\ddj=zeros(NDOF); 

dxdj=zeros(l.NDOF)'; 

dxddj=zeros(l.NDOF)'; 

forn=l:NNODES 

ny=4*(n-l)+l; 

nz=4*(n-l)+2; 

dxdj(ny)=0.01; 

dxdj(nz)=0.01; 

dxddj(ny)=0.01; 

dxddj(nz)=0.01; 

dFdxdj  l(ny)=... 

...(-0.5*rho*By(n)*Cdly(n)*(xdj(ny)+dxdj(ny))*abs(xdj(nv)+dxdj(ny))*. 

...Lh(n)- 
0.5*rho*Dapy(n)*Cd2y(n)*(xdj(ny)+dxdj(ny))*abs(xdj(ny)+dxdj(ny))*... 

...Lapy(n)-Fvel(ny))/dxdj(ny); 

dFdxdj  l(nz)=... 

...(-5*rho*Bz(n)*Cdlz(n)*(xdj(nz)+dxdj(nz))*abs(xdj(nz)+dxdj(nz))*... 

...Lh(n)- 
0.5*rho*Dapz(n)*Cd2z(n)*(xdj(nz)+dxdj(nz))*abs(\dj(nz)+dxdj(nz))*... 

..  Lapz(n)-Fvel(nz))/dxdj(nz); 

dFdxddj  l(ny)=(-rho*(pi*By(n)A2/4)*Caly(n)*Lh(n)*... 

...(xddj(ny)+dxddj(ny))-rho*Aapy(n)*Ca2y(n)*Lapy(n)*... 

...(xddj(ny)+dxddj(ny))-Facc(ny))/dxddj(ny); 

dFdxddj  l(nz)=(-rho*(pi*Bz(n)A2/4)*Calz(n)*Lh(n)*... 

...(xddj(nz)+dxddj(nz))-rho*Aapz(n)*Ca2z(n)*Lapz(n)*... 

...(xddj(nz)+dxddj(nz))-Facc(nz))/dxddj(nz); 
end 
for  n=l:  NDOF 

dFdxdj(n.n)-dFdxdj  l(n); 

dFdxddj  ( n.  n  )=dFdxddj  1  ( n ) ; 
end 
% 

%  Calculate  residuals 

flj=M*(cl*(xj-xt)-c2*xdt-c3*xddt)+C*(c4*(xj-xt)-c5*xdt+c6*xddt*dt)+... 
...K*(xj-xt)-Floadj+Floadt; 
f2j=xdj-c4*(xj-xt)+(c5-l)*xdt-c6*dt*xddt; 
f3j=xddj-cl*(xj-xt)+c2*xdt+c7*xddt; 
% 

%  Solve  for  the  displacement  correction  (dx) 
% 

Al=(cl*M+c4*C+K-cl*dFdxddj-c4*dFdxdj-dFdxj): 
Alinv-inv(Al). 

A2=-flj-dFdxdj*f2j-dFdxddj*Bj; 
dx=(Alinv*A2); 
% 

%  Calculate  the  next  iteration's  displacements,  velocities,  and  accelerations 
% 

xj=xj+dx; 
xdj=xdj+(c4*dx)-f2j; 
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\ddj=xddj+(cl*d\)-Oj; 

% 

%  Check  to  see  if  all  Norms  of  residuals  are  less  than  tolerance  for 

%  the  current  iteration 

%  (If  so.  go  to  the  next  time  step) 

Normflj=norm(flj); 

Normf2j=norm(f2j); 

Normf3j=norm(f3j); 

if  Normflj<=tol  &  Normf2j<=tol  &  Normf3j<=tol 

/o 

%  Save  current  time,  displacements,  velocities. 

%  and  accelerations  (at  specified  nodes)  for  later  plotting 

timefi)^: 

%  DisplacementsA  elocities/accelerations 

for  n=  1:4 

Dh(n)=(N(n)-l)*4+2; 

end 

Xl(i)=xj(Dh(l)); 

Xdl(i)=xdj(Dh(l)): 

Xddl(i)=xddj(Dh(l)); 

X2(i)=xj(Dh(2)): 

Xd2(i)-xdj(Dh(2)): 

Xdd2(i)=xddj(Dh(2)); 

X3(i)=xj(Dh(3)); 

Xd3(i)=xdj(Dh(3)): 

Xdd3(i)=xddj(Dh(3)): 

X4(i)=xj(Dh(4)): 

Xd4(i)=xdj(Dh(4)): 

Xdd4(i)=xddj(Dh(4)); 

% 

Dhl9=(19-l)*4+2: 

Dhl8=(18-l)*4+2: 

Fvell8(i)=Fvel(Dhl8): 

Faccl8(i)=Facc(Dhl8); 

Floadl8(i)-Floadj(Dhl8): 

Fvell9(i)=Fvel(Dhl9); 

Faccl9(i)=Facc(Dhl9); 

Floadl9(i)=Floadj(Dhl9): 

% 

%  Save  the  current  time  step  displacements,  velocities,  and 

%  accelerations,  and  forces  for  use  in  the  next  time  step 

xt=xj; 

\dt=xdj: 

xddt=xddj: 

Floadt=Floadj; 

break,  end 
end 
if  j==jlimit.  break,  end 
end 
if  j==jlimit.  break,  end 
if  t>=tf.  break,  end 
end 
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%  Subroutine  BUBBLE. M 

% 

%  This  subroutine  performs  the  following: 

%  1.  Computes  the  time  histories  of  the  explosion  bubble  radius  (rad) 

%  and  bubble  depth  (depth)  by  integration  of  the  differential 

%  equations  governing  bubble  dynamics. 

%  Computation  is  earned  out  for  specified  bubble  initial  conditions 

%  (i.e.  initial  depth  of  the  charge  (Dchg)  and  weight  of  the 

%  charge  (Wchg)). 

% 

% 

function[TM.BR.BD.BRD.BDD.BRDD.BDDD.DEL.btime.brad.bdpth.bvel.Ll.Tl|= 

. .  bubble(Dchg.Wchg.tmax) 

% 

%  Define  constants 

% 

%  Drag  coeff.  for  bubble 

Cd=2.25; 

%  Initial  bubble  head  (ft) 

DO=Dchg+33; 

%  Adiabatic  gas  exponent  for  TNT 

gl=1.25; 

%  Constant  for  TNT 

kl=(0.0552)*D()A0.25; 

%  Length  scale  constant 

Ll=13.67*(Wchg/D())A(l/3); 

%  Time  scale  constant 

Tl=2.94*(WchgA(l/3)/D0A(5/6)); 

%  Misc.  constants 

cl=33/Ll; 

c2=D0/Ll; 

c3=(gl-l)*kl; 

c4=3*gl+l. 

% 

%  Define  initial  and  final  conditions 

% 

%  Non-dim.  bubble  radius 

xOklA(4/3)*(l+(4*klA4)/3): 

%  Non-dim.  bubble  depth 

y0=c2; 

%  Rate  of  change  of  bubble  radius 

xdO=0; 

%  Rate  of  change  of  bubble  depth 

ydO=0; 

%  Non-dim.  start  time  (i.e.  0  seconds) 

TO-0.0; 

%  Non-dim.  end  time  (i.e.  tmax  seconds) 

Tf=tmax/Tl; 

%  Non-dimensional  time  step 

dT=(5.0*xOA4)/Tf; 

%  Integration  time  step  parameter 

tsp=10.; 

% 

%  Perform  Runge-Kutta  time  iteration 
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"/o 

T=TO; 

x=xO; 

y%0; 

xd=xdO: 

yd=ydO; 

% 

M=3000;  %  defined  maximum  number  of  time  steps  (for  control) 

% 

for  m=l:M 

%  Call  the  differential  eqn.  for  bubble  radius 

[x.xd]=radfn(dT.x.xd.\  .yd.c  1  .c2.c3.c4.Cd): 

xdd=-1.5/(l-x/(2*(y-ci)))*((xdA2/x)*(l-2*x/(3*(y-cl)))-ydA2/(6*x)+y/(x*c2)-c3/(xAc4)+. 

. .  .(x/(4*(y-c  1  )A2))*(.\d*yd/3  -x/c2+Cd*ydA2/4)); 

%  Call  the  differential  eqn.  for  bubble  depth 

Lv.yd|=dpthfn(dT.x.xd.y.yd.xdd.cl.c2.c3.c4.Cd); 

ydd=-3*(l/c2+xd*vd/x-(Cd*vdA2)/(4*x)+(x/(4*(y<l)A2))*(3*xdA2+x*xdd)); 

% 

%  Save  the  current  bubble  characteristics  for  plotting  and  later  use 

% 

TM(m)=T;  %  non-dim.  time 

BR(m)=x;  %  non-dim  bubble  radius 

BD(m)=y;  %  non-dim.  bubble  depth 

BRD(m)=xd;  %  non-dim.  rate  of  change  of  bubble  radius 

BDD(m)=yd;  %  non-dim.  rate  of  change  of  bubble  depth  (bubble  velocity) 

BRDD(m)=.\dd; 

BDDD(m)=ydd, 

DEL(m)=(v-cl);  %  non-dim.  suface  pressure 

% 

btime(m)=T*Tl;  %  dimensional  time 

brad(m)=x*Ll;  %  dimensional  bubble  radius 

bdpth(m)=(y-cl)*Ll;  %  dimensional  bubble  depth 

bvel(m)=yd*Ll/Tl;  %  dimensional  bubble  velocity 

% 

%  Go  to  next  time  increment 

x3=xA3; 

x3d=3.0*(xA2)*xd; 

x3dd=3.0*x*(x*xdd+2.0*xdA2); 

dT=2.0/(abs(.\3dd)*3.0*tsp); 

ifdT<0.001.dT=0.0()l; 
elseif  dT>0.025.  dT=0.025: 
end 

T=T+dT; 

if  T>Tf,  break 
end 
end 
% 

%  Plot  the  results 
% 

figure(l) 

subplot(3.1.1).  plot(btime.brad) 
title( 'Bubble  radius  vs.  time') 
xlabelfTime  (sec)') 
ylabeK'Bubble  radius  (ft)') 
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subplot(3. 1.2).  plotfbtime.-bdpth) 
titleCBubble  depth  vs.  time') 
xlabeK'Time  (sec)') 
ylabeK 'Bubble  depth  (-  ft)') 
subplot(3. 1.3).  plot(btime.-bvel) 
title( 'Upward  bubble  velocity  vs.  time') 
xlabeK'Time  (sec)') 
ylabeK 'Upward  bubble  velocity  (ft/s)') 
set(gcf.'PaperPosition'.|1.0.  2.0.  6.5.  8.5]) 
%  print 
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%  RADFN.M 

% 

%  Function  M-file  for  calculation  of  the  bubble  radius 

%  (i.e.  contains  the  bubble  radius  governing  equation  for  use  in  BUBBLE  M) 

% 

function  [x.xdj=radfn(dT.x.xd.y.vd.c  1  .c2.c3.c4.Cd) 

akl=dT*rfunc(x.xd.y.yd.cl.c2.c3.c4.Cd); 

ak2=dT*rfunc(x+dT*xd/2.0.xd+akl/2.0.y.yd,cl.c2.c3.c4.Cd); 

ak3=dT*rfunc(x+dT*xd/2.0+dT*akl/4  0.xd+ak2/2.().y.yd.cl.c2.c3.c4.Cd); 

ak4=dT*rfunc(x+dT*xd+dT*ak2/2.0.xd^ak3.y.yd.cl.c2.c3.c4.Cd); 

x=x+dT*xd+dT*(ak  1  +ak2+ak3  )/6.0; 

xd=xd+(akl+2.()*ak2+2.0*ak3+ak4)/6.0; 


%  RFUNC.M 

% 

function  ak=rfunc(x,xd,y,yd.cl,c2.c3.c4.Cd) 

ak=-1.5/(l-x/(2*(y<l)))*((xdA2/x)*(l-2*x/(3*(y-cl)))-ydA2/(6*x)+y/(x*c2)-c3/(xAc4)+(x/(4*(y- 

cl)A2))*(xd*yd/3-x/c2+Cd*ydA2/4)); 
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%  DPTHFN.M 

% 

%  Function  M-file  for  calculation  of  the  bubble  depth 

%  (i.e.  contains  the  bubble  depth  governing  equation  for  use  in  BUBBLE. M) 

/o 

function  |v.vd]=dpthfn(dT.x.xd.y.vd.xdd.cl.c2.c3.c4.Cd) 

akl=dT*dfunc(x.xd.y.vd.\dd.cl.c2.c3.c4.Cd); 

ak2-dT*dfunc(.\.xd.y+dT*yd/2.0.vd+akl/2.0.xdd.cl.c2.c3.c4.Cd): 

ak3=dT*dfunc(x.xd.y+dT*yd/2.0+dT*akl/4.().yd+ak2/2.0.xdd.cl.c2.c3.c4.Cd); 

ak4=dT*dfunc(x.xd.y+dT*yd+dT*ak2/2.0.yd+ak3.xdd.cl.c2.c3.c4.Cd): 

y=y+dT*yd+dT*(akl+ak2+ak3)/6.0; 

yd=yd+(akl+2.0*ak2+2.0*ak3+ak4)/6.0; 


%  DFUNC.M 

% 

% 

% 

function  ak=dfunc(x.xd.y.yd.xdd.c  1  .c2.c3.c4.Cd) 

ak=-3*(l/c2+xd*yd/x-(Cd*ydA2)/(4*x)+(x/(4*(y-cl)A2))*(3*xdA2+x*xdd)); 
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%  BUB  WHIP  M 

% 

% 

%  This  subroutine  performs  the  following: 

%  1.  Calculates  the  response  (nodal  point  displacements,  velocities,  accelerations) 

%  of  the  submarine  during  a  vibration  with  initial  displacements  of  zero 

%  and  subjected  to  fluid  loading  produced  by  a  pulsating  explosion  bubble. 

%  Response  is  calculated  using  the  Newmark  time  integration  algorithm. 

%  with  iteration  conducted  at  each  time  step  using  a  Newton-Raphson  method 

%  2.  Calculates  the  local  axial  strains  at  specified  hull  locations  using  simple 

%  beam  theorv  (used  for  comparison  with  model  test  data). 

% 

% 

%  Utilize  the  following: 

%  -K.C.M  calculated  in  HULL.M 

%  -Hydrodynamic  coefficients  provided  in  IN  M 

%  -Parameters  for  Newmark  algorithm  provided  in  IN  M 

%  -Bubble  characteristics  (radius,  depth,  etc  )  calculated  in  BUBBLE  M 

% 

%  Define  nodes  to  use  for  response  calculations 

N(l)=NODEl; 

N(2)=NODE2; 

N(3)=NODE3; 

N(4)=NODE4; 

%  Define  hull  elements  to  use  for  axial  strain  calculations 

E11=ELEMENT1; 

E12=ELEMENT2; 

E13=ELEMENT3; 

% 

%  Define  the  horizontal  (hull)  degree  of  freedom  where  internal  structure/absorber  acts 

nnabs=(  nabs- 1  )*4+2 ; 

% 

%  Define  parameters  for  internal  structure/absorber 

ka=ma*(ww(8)A2); 

mabs=zeros(  1  .NDOF)'; 

kabs-zeros(l.NDOF)'; 

mabs(nnabs)=ma: 

kabs(nnabs)=ka; 

% 

%  Set  the  initial  and  final  times  (sec)  for  Newmark  integration,  and  limit 

%  on  the  number  of  Newton-Raphson  iterations 

t=0; 

tf=tmax, 

jlimit=40; 

% 

%  Initialize  the  loading  forces  on  the  hull 

%  Loading  forces  related  to  displacement 

Fdisp=zeros(l.NDOF)*: 
%  Loading  forces  related  to  velocity 

Fvel=zeros(l.NDOF)'; 
%  Loading  forces  related  to  acceleration 

Facc=zeros(l.NDOF)'; 
%  Calculate  (constant)  static  forces  (i.e.  hydrostatic/buoyancy  and  weight) 
Fgrav=zeros(l.NDOF)'; 
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Fhs-zeros(l.NDOF)'; 
forn=l:NNODES 

ny=(n-l)*4+l; 

%  Hydrostatic/buoyancy  force  (lbf) 

Fhs(ny)=bh(n)*2240: 

%  Static  force  of  gravity/weight  (lbf) 

Fgrav(ny)=-vvh(n)*2240; 
end 

%  Total  initial  loading  forces 
Floadt=Fdisp+Fvel+Facc+Fhs+Fgrav; 
% 

%  Define  initial  conditions  (displacements/velocities/accelerations) 
%  Hull 

xt=zeros(l.NDOF)'; 
xdt=zeros(l.NDOF)'. 
xddt=zeros(l.NDOF)': 
%  Internal  absorber 
xat=0; 
xadt-0; 
xaddt=0; 
% 

%  Perform  the  Nevvmark  time  integration 
% 

%  Step  through  time 
% 

%  Define  constants  to  be  used  in  Newton-Raphson  iterations 
cl=(l/(betap*dtA2)): 
c2=(l/(betap*dt)); 
c3=(l/(2*betap)); 
c4=(gammap/(betap*dt)); 
c5=gammap/betap; 
c6=(  1  -gammap/(  2  *betap ) ) ; 
c7=(l/(2*betap)-l); 
% 

dx=zeros(l.NDOF)'; 
% 
for  i=l:199; 

t=t+dt 

% 

%  Set  initial  displacments.velocities.accelerations.  and  loading  forces  for 

%  the  current  time  step  (equal  to  the  final  for  the  last  time  step) 

xj=xt; 

xdj=xdt; 

xddj=xddt; 

Floadj=Floadt; 

xaj=xat; 

xadj=xadt; 

xadds=xaddt; 

% 

%  Perform  Newton-Raphson  iterations  at  current  time  step 

for  j=l  jlimit 
% 

\addj=xadds; 
%  Calculate  current  iteration  absorber  displacement/velocity 
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xaj=xat+xadt*dt+xaddt*(dtA2/2)+betap*(xaddj-xaddt)*dtA2; 

xadj=xadt+xaddt*dt+gammap*(xaddj-xaddt)*dt; 

% 

%  Calculate  the  loading  forces  for  the  current  iteration 

% 

%  Call  function  "fluid"  to  provide  free  field  fluid  velocities  and 

%  accelerations  at  each  node  (at  the  current  iteration  node  location) 

[uy.uz.uyd.uzdl=fluid(xj.xchg.ychg.zchg.Dship.t.TM.BR.BD.BRD.BDD.BRDD.BDDD. 
...DEL.Ll.Tl.btime.bdpth.NDOF.NNODES.Lh); 
% 
forn=l:NNODES 

ny=4*(n-l)+l; 

nz=4*(n-l)+2; 
%  Loading  forces  related  to  displacement 

Fdisp(nz)=kabs(nz)*(xaj-xj(nz)); 
%  Loading  forces  related  to  velocity 

Fvel(ny)=0.5*rho*By(n)*Cdly(n)*(uy(n)-xdj(ny))*... 

...abs(uy(n)-\dj(ny))*Lh(n)+()  5*rho*Dapy(n)*Cd2y(n)*(uy(n)- 

...xdj(ny))*abs(uy(n)-\dj(ny))*Lapy(n); 

Fvel(nz)=0.5*rho*Bz(n)*Cdlz(n)*(uz(n)-xdj(nz))*... 

...abs(uz(n)-xdj(nz))*Lh(n)+()5*rho*Dapz(n)*Cd2z(n)*(uz(n)- 

...xdj(nz))*abs(uz(n)-xdj(n/.))*Lapz(n); 
%  Loading  forces  related  to  acceleration 

Facc(ny)=rho*(pi*By(n)A2/4)*Cmly(n)*Lh(n)*(uyd(n)- 

...xddj(ny))+rho*Aapy(n)*Cm2y(n)*Lapy(n)*(uyd(n)- 
...xddj(ny))+rho*(pi*By(n)/N2/4)*Lh(n)*xddj(ny)+rho*Aapy(n)*Lapy(n)*xddj(ny); 

Facc(nz)=rho*(pi*Bz(n)A2/4)*Cmlz(n)*Lh(n)*(uzd(n)- 

...xddj(nz))+rho*Aapz(n)*Cm2z(n)*Lapz(n)*(uzd(n)- 
...xddj(nz))+rho*(pi*Bz(n)A2/4)*Lh(n)*xddj(nz)+rho*Aapz(n)*Lapz(n)*xddj(nz); 
end 
% 

%  Total  loading  force  for  the  current  iteration 
Floadj=Fdisp+Fvel+Facc+Fhs+Fgrav; 
% 

%  Calculate  Jacobian  matrices  for  each  of  the  loading  forces 
%  (square  diagonal  matrices  for  this  formulation) 
% 

dFdxj  1  =zeros(  1  ,NDOF)'; 
dFdxdj  1  =zeros(  1  .NDOF)'; 
dFdxddj  1  =zeros(  1  .NDOF)'; 
dFdxj=zeros(NDOF); 
dFdxdj=zeros(NDOF); 
dFdxddj=zeros(  NDOF); 
dxj=zeros(l.NDOF)'; 
dxdj=zeros(l,NDOF)'; 
dxddj=zeros(  LNDOF)'; 
forn=l:NNODES 

ny=4*(n-l)+l; 

nz=4*(n-l)+2; 

dxj(nz)=0.01; 

dxdj(ny)=0.01; 

dxdj(nz)=0.01; 

dxddj(ny)=0.01; 
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dxddj(nz)=0.01; 

dFdxjl(nz)=(kabs(nz)*(xaj-(xj(nz)+dxj(nz)))-Fdisp(nz))/dxj(nz): 
dFd\djl(ny)=(()5*rho*By(n)*Cdly(n)*(uy(n)-(.\dj(ny)+d\dj(ny)))*abs(uy(n)- 
...(xdj(ny)+d.\dj(ny)))*Lh(n)+0.5*rho*Dapy(n)*Cd2y(n)*(uy(n)- 
...(xdj(ny)+dxdj(ny)))*abs(uy(nH\dj(ny)+dxdj(ny)))*Lapy(n)- 
...Fvel(ny))/dxdj(ny); 

dFd\djl(nz)=(0  5*rho*Bz(n)*Cdlz(n)*(uz(n)-(xaj(nz)+d\dj(nz)))*abs(uz(n)- 
...(xdj(nz)+dxdj(nz)))*Lh(n)+0.5*rho*Dapz(n)*Cd2z(n)*(uz(n)- 
...(xdj(nz)+dxdj(nz)))*abs(uz(n)-(xdj(nz)+dxdj(nz)))*Lapz(n)- 
...Fvel(nz))/dxdj(nz); 

dFd\ddjl(ny)=(rho*(pi*By(n)A2/4)*Cmly(n)*Lh(n)*(uyd(n)- 
...(xddj(ny)+dxddj(ny)))+rho*Aapy(n)*Cm2y(n)*Lapy(n)*(uyd(n)- 
(xddj(ny)+dxddj(ny)))+rho*(pi*By(n)A2/4)*Lh(n)*(\ddj(ny)+d\ddj(ny))+rho*Aapy(n)* ... 
...Lapy(n)*(xddj(ny)+dxddj(ny))-Facc(ny))/d\ddj(ny); 
dFdxddjl(nz)=(rho*(pi*Bz(n)A2/4)*Cmlz(n)*Lh(n)*(uzd(n)- 
...(xddj(nz)+dxddj(nz)))+rho*Aapz(n)*Cm2z(n)*Lapz(n)*(uzd(n)- 
,(xddj(nz)+dxddj(nz)))+rho*(pi*Bz(n)A2/4)*Lh(n)*(.\ddj(nz)+dxddj(nz))+rho*Aapz(n)*... 
..Lapz(n)*(xddj(nz)+dxddj(nz))-Facc(nz))/dxddj(nz); 

end 

forn=l:NDOF 

dFdxj(n.n)=dFdxjl(n); 
dFdxdj(n.n)=dFdxdj  l(n); 
dFdxddj(n.n)=dFd\ddj  l(n). 

end 

% 

%  Calculate  residuals 

nj=M*(cl*(xj-xt)-c2*xdt-c3*xddt)+C*(c4*(xj-xt)-c5*xdt+c6*xddt*dt)+... 

. .  .K*(xj-xt)-Floadj+Floadt; 

f2j=xdj-c4*(xj-xt)+(c5-l)*xdt-c6*dt*xddt; 

Oj=xddj-c  1  *(xj-xt)+c2*xdt+c7*xddt: 

% 

%  Solve  for  the  displacement  correction  (dx) 

% 

Al=(c  1  *M+c4*C+K-c  1  *dFdxddj-c4*dFdxdj-dFdxj); 

Alinv=inv(Al); 

A2=-flj-dFdxdj*f2j-dFdxddj*f3j; 

dx=(Alinv*A2); 

% 

%  Calculate  the  next  iteration's  displacements,  velocities,  and  accelerations 

% 

xj=xj+dx; 

xdj=xdj+(c4*dx)-f2j; 

xddj=xddj+(cl*dx)-f3j; 

% 

%  Required  acceleration  of  absorber 

if  any(mabs)~=0.  xadds=-Fdisp(nnabs)/mabs(nnabs); 

else.  xadds=0; 

end 

% 

%  Check  to  see  if  all  Norms  of  residuals  are  less  than  tolerance  for 

%  the  current  iteration 

%  (If  so.  go  to  the  next  time  step) 

Normflj=norm(flj); 

Normf2j=norm(f2j); 
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Normf3j=norm(f3j); 
Normabs=norm(  xadds-xaddj ); 

if  Normflj<=tol  &  Normf2j<=tol  &  Normf3j<=tol  &  Normabs<=tol 
% 

%  Calculate  axial  stress  at  specified  hull  stations 
%  (using  simple  beam  theory) 
% 

qy=zeros(l.NE)': 
qz=zeros(l.NE)'; 
Vy=zeros(l.NE)'; 
Vz=zeros(l.NE)'; 
Myl=zeros(l,NE)': 
Mzl=zeros(l.NE)*: 
My=zeros(l.NE)'; 
Mz=zeros(I.NE)'; 
fore=l:NE 

%  Net  vertical  or  horizontal  force  at  element  e 

ey=(e-l)*4+l; 

ez=(e-l)*4+2; 

qy(e)=Floadj(ey)-M(ey.ey)*xddj(ey): 

qz(e)=Floadj(ez)-M(ez.ez)*.\ddj(ez): 

% 

%  Vertical  or  horizontal  shear  force  at  element  e 

Vy(e)=sum(qy); 

Vz(e)=sum(qz); 

% 

%  Vertical  or  horizontal  bending  moment  at  element  e 

Myl(e)=Vy(e)*Le(e); 

Mzl(e)=Vz(e)*Le(e); 

My(e)=sum(Myl);     %  BM  about  z-axis 

Mz(e)=sum( Mz  1 );  %  BM  about  y-axis 

% 

%  Axial  stress  at  crown  of  element  e  (positive  is  tension) 

sigc(e)=-My(e)*Yc(e)/laz(e); 

%  Axial  stress  at  keel  of  element  e 

sigk(e)=My(e)*Yk(e)/Iaz(e); 

%  Axial  stress  at  port-side  fiber  (port=pos  z)  (positive  is  tension) 

sigp(e)=-Mz(e)*Zp(e)/Iay(e); 

%  Axial  stress  at  stbd-side  fiber 

sigs(e)=Mz(e)*Zs(e)/Iay(e); 

% 

%  Axial  strains 

EPSc(e)=sigc(e)/E; 

EPSk(e)=sigk(e)/E; 

EPSp(e)=sigp(e)/E; 

EPSs(e)=sigs(e)/E; 

% 
end 

%  Save  current  time,  displacements  and  v  elocities 
%  (at  specified  nodes)  and  axial  strains 
%  (at  specified  elements)  for  later  plotting...  also  save 
%  current  horizontal  fluid  velocity  and  acceleration,  and 
%  (at  nodes  10  and  20)  and  all  calculated  horizontal  forces 
%  (at  nodes  10  and  20)  for  later  plotting. 
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time(i)=t; 

%  Vertical/horizontal  displacements/velocities 

%  at  specified  nodes 

for  n=  1:4 

Dv(n)=(N(n)-l)*4+l; 
Dh(n)=(N(n)-l)*4+2; 

end 

xvl(i)=xj(Dv(l)); 

xhl(i)=xj(Dh(I)); 

xdvl(i)=xdj(Dv(l)); 

xdhl(i)=xdj(Dh(l)); 

xv2(i)=xj(Dv(2)); 

xh2(i)=xj(Dh(2)); 

xdv2(i)=xdj(Dv(2)); 

xdh2(i)=xdj(Dh(2)); 

xv3(i)=xj(Dv(3)); 

xh3(i)=xj(Dh(3)); 

xdv3(i)=xdj(Dv(3)); 

xdh3(i)=xdj(Dh(3)); 

xv4(i)=xj(Dv(4)); 

xh4(i)=xj(Dh(4)); 

xdv4(i)=xdj(Dv(4)); 

xdh4(i)=xdj(Dh(4)); 

% 

xabs(i)=xaj; 

xdabs(i)=xadj; 

xddabs(i)=xaddj; 

% 

%  Horizontal  hull/fluid  motions  and 

%  external  forces 

xddh2(i)=xddj(Dh(2)); 

xddh3(i)=xddj(Dh(3)); 

uh3(i)=uz(N(3)); 

udh3(i)=uzd(N(3)); 

Fvelh3(i)=Fvel(Dh(3)); 

Facch3(i)=Facc(Dh(3)); 

Floadh3(i)=Floadj(Dh(3)); 

% 

xddh4(i)=xddj(Dh(4)); 

uh4(i)=uz(N(4)); 

udh4(i)=uzd(N(4)); 

Fvelh4(i)=Fvel(Dh(4)); 

Facch4(i)=Facc(Dh(4)); 

Floadh4(i)=Floadj(Dh(4)); 

% 

%  Axial  strains  (keel,  crown,  port,  stbd)  at 

%  at  specified  elements 

strainkl(i)=EPSk(Ell); 

straincl(i)=EPSc(Ell); 

strainpl(i)=EPSp(Ell); 

strains  l(i)=EPSs(Ell); 

straink2(i)=EPSk(E12); 

strainc2(i)=EPSc(E12); 

strainp2(i)=EPSp(E12); 
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strains2(i)=EPSs(E12); 

straink3(i)=EPSk(E13); 

strainc3(i)=EPSc(E13); 

strainp3(i)=EPSp(E13); 

strains3(i)=EPSs(E13); 

% 

%  Save  the  current  time  step  displacements,  velocities,  and 

%  accelerations,  and  forces  for  use  in  the  next  time  step 

xt=xj; 

xdt=xdj: 

xddt=xddj: 

Floadt=Floadj; 

xat=xaj; 

xadt=xadj; 

\addt=xaddj; 

break,  end 
end 
if  j==jlimit.  break,  end 
end 
if  j==jlimit.  break,  end 
if  t>=tf.  break,  end 
end 
% 
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%  FLUID  M 

% 

%  Calculates  the  free  field  fluid  velocities  and  accelerations  (at  each  node) 

% 

function[uy.uz.uyd.uzd]=fluid\el(xj.xchg.vchg.zchg.Dship.t.TM.BR.BD.BRD.BDD 

...BRDD.BDDD.DEL.Ll.Tl.btime.bdpth.NDOF.NNODES.Lh) 

% 

%  Non-dimensional  bubble  parameters 

TI=t/Tl; 

a=interpl(TM.BR.TI); 

s=interpl(TM.BRD.TI); 

sd=interpl(TM.BRDD.TI); 

l=interpl(TM.BDD.TI). 

ld=interpl(TM.BDDD.TI); 

d=interpl(TM.DEL.TI): 

%  Source  strengths 

el=(LlA3)*(aA2)*s/Tl; 

e2=-(LlA4/(2*Tl))*(aA3*l); 

eld=(LlA3/TlA2)*(aA2*sd+2*a*sA2): 

e2d=-(LlA4/(2*TlA2))*(aA3*ld+3*aA2*l*s); 

%  Vertical  bubble  velocity 

vm=-l*Ll/Tl; 

%  Bubble  depth 

ti=t; 

bdepth=interpl(btime.bdpth.ti): 

% 

forn-l:NNODES 

ny=(n-l)*4+l; 

nz=(n-l)*4+2; 

xnode=(n-l)*Lh(n); 

ynode=xj(ny): 

znode=xj(nz); 

%  vr  and  hr  (vertical  and  horiz  relative  distance)  and  angle  of  flow 

vr=bdepth-Dsrup+ynode; 

hr=((xnode-xchg)A2+(znode-zchg)A2)A0.5; 

theta=atan((xnode-xchg)/(znode-zchg)); 

% 

%  Free-field  fluid  velocities/accelerations  (vertical  and  horizontal) 

%  Horizontal  and  vertical  free-field  fluid  velocities  and  accelerations 

dis=(vrA2+hrA2)A0.5; 

uh=(el*hr/disA3)+(3*e2*vr*hr/disA5); 

uv=(el*vr/disA3)+(3*e2*vrA2/disA5)-(e2/disA3); 

uhd=(eld*hr/disA3)+(3*e2d*hr*vr/disA5)+... 

...vm*(3*el*hr*vr/disA5-(3*e2*hr/disA5)*(l-(5*vrA2/disA2))); 

uvd=(eld*vr/disA3+3*e2d*vrA2/disA5-e2d/disA3)-... 

...vm*(el/disA3-3*el*vrA2/disA5+9*e2*vr/disA5-15*e2*vrA3/disA7); 

% 

%  Y  and  Z  directed  free-field  fluid  velocities  and  accelerations 

uy(n)=uv; 

uz(n)=uh*cos(theta); 

uyd(n)=uvd; 

uzd(n)=uhd*cos(theta); 
end 
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APPENDIX  B 


•DRY"  WHIPPING  RESPONSE  AT  SPECIFIED  NODES 
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Figure  B-l.  Horizontal  response  at  node  7,  given  initial  displacments  equal  to  a  scaled 
horizontal  mode  shape  vector  for  mode  1  flexure 
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Horizontal  displacement  vs.  time  for  node  10 
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Figure  B-2.  Horizontal  response  at  node  10,  given  initial  displacments  equal  to  a  scaled 
horizontal  mode  shape  vector  for  mode  1  flexure 
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Figure  B-3.  Horizontal  response  at  node  20,  given  initial  displacments  equal  to  a  scaled 
horizontal  mode  shape  vector  for  mode  1  flexure. 
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Figure  B-4.  Horizontal  response  at  node  7,  given  initial  displacments  equal  to  a  scaled 
horizontal  mode  shape  vector  for  mode  2  flexure 
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Figure  B-5.  Horizontal  response  at  node  10,  given  initial  displacments  equal  to  a  scaled 
horizontal  mode  shape  vector  for  mode  2  flexure. 
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Figure  B-6    Horizontal  response  at  node  20,  given  initial  displacments  equal  to  a  scaled 
horizontal  mode  shape  vector  for  mode  2  flexure. 
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APPENDIX  C 


"WET"  WHIPPING  RESPONSE  AT  SPECIFIED  NODES 
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Figure  C-l.  Horizontal  response  at  node  7,  given  initial  displacments  equal  to  a  scaled 
horizontal  mode  shape  vector  for  mode  1  flexure. 
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Figure  C-2.  Horizontal  response  at  node  10,  given  initial  displacements  equal  to  a  scaled 
horizontal  mode  shape  vector  for  mode  1  flexure. 
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Figure  C-3.  Horizontal  response  at  node  20,  given  initial  displacments  equal  to  a  scaled 
horizontal  mode  shape  vector  for  mode  1  flexure. 
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Figure  C-4.  Horizontal  response  at  node  7,  given  initial  displacments  equal  to  a  scaled 
horizontal  mode  shape  vector  for  mode  2  flexure. 
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Figure  C-5    Horizontal  response  at  node  10,  given  initial  displacments  equal  to  a  scaled 
horizontal  mode  shape  vector  for  mode  2  flexure. 
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Horizontal  displacement  vs.  time  for  node  20 
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Figure  C-6    Horizontal  response  at  node  20,  given  initial  displacments  equal  to  a  scaled 
horizontal  mode  shape  vector  for  mode  2  flexure. 
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APPENDIX  D 


WHIPPING  RESPONSE  FOR  HULL  SUBJECTED  TO  CHARGE  DESIGNED  TO 
EXCITE  HORIZONTAL  MODES  1   AND  2  FLEXURAL  WHIPPING 
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Figure  D- 1    Horizontal  response  at  node  7,  for  hull  subject  to  bubble  pulse  designed  to 
excite  horizontal  mode  1  flexural  whipping. 
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Figure  D-2    Horizontal  response  at  node  10,  for  hull  subject  to  bubble  pulse  designed  to 

excite  horizontal  mode  1  flexural  whipping. 


154 


Horizontal  displacement  vs.  time  for  node  20 


Time  (sec) 


o 

c/) 


>>   o 


o 
o 

a) 

> 


Horizontal  velocity  vs.  time  for  node  20 


1  1.5 

Time  (sec) 


2.5 


-1000 


Horizontal  acceleration  vs.  time  for  node  20 


0.5 


1  1.5 

Time  (sec) 


Figure  D-3.  Horizontal  response  at  node  20,  for  hull  subject  to  bubble  pulse  designed  to 

excite  horizontal  mode  1  flexural  whipping. 
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Figure  D-4    Horizontal  response  at  node  7,  for  hull  subject  to  bubble  pulse  designed  to 
excite  horizontal  mode  2  flexural  whipping. 
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Figure  D-5.  Horizontal  response  at  node  10,  for  hull  subject  to  bubble  pulse  designed  to 

excite  horizontal  mode  2  flexural  whipping. 
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Figure  D-6.  Horizontal  response  at  node  20,  for  hull  subject  to  bubble  pulse  designed  to 

excite  horizontal  mode  2  flexural  whipping. 
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APPENDIX  E 


WHIPPING  RESPONSE  FOR  HULL  SUBJECTED  TO  CHARGE  DESIGNED  TO 
EXCITE  HORIZONTAL  MODE  2  FLEXURAL  WHIPPING  AND  INCLUDING 

INTERNAL  MASS/ ABSORBER 
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Horizontal  displacement  vs.  time  for  node  7 
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Figure  E-l    Horizontal  response  at  node  7,  for  hull  subject  to  bubble  pulse  designed  to 
excite  honzontal  mode  flexural  whipping  and  including  200  LT  "tuned"  mass  damper  at 

node  7. 
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Horizontal  displacement  vs.  time  for  node  10 
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Figure  E-2.  Horizontal  response  at  node  10,  for  hull  subject  to  bubble  pulse  designed  to 
excite  horizontal  mode  flexural  whipping  and  including  200  LT  "tuned"  mass  damper  at 

node  7. 
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Figure  E-3.  Horizontal  response  at  node  20,  for  hull  subject  to  bubble  pulse  designed  to 
excite  horizontal  mode  flexural  whipping  and  including  200  LT  "tuned"  mass  damper  at 

node  7. 


162 


Horizontal  displacement  vs.  time  for  internal  absorber 
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Figure  E-4    Horizontal  response  of  200  LT  "tuned"  mass  located  at  node  7 
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APPENDIX  F 


NATURAL  SLOSHING  FREQUENCIES  AND  SLOSHING  MODAL  MASS  FOR 
LATERAL  SLOSHING  OF  RECTANGULAR  TANKS 
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Figure  F-l .  Natural  sloshing  frequency  and  sloshing  modal  mass  for  lateral  sloshing  of  a 

(10'h  x  10*w  x  lO'l)  rectangular  tank. 


165 


o20 

03 

o 

§  10 

cr 
d) 

l_ 

z    0 


0 


2  3 

Sloshing  mode  (n) 


c 
o 

CD 

£ 

|  40 

o 

E 

°20 

.c 

CD 


2  3  4 

Sloshing  mode  (n) 


Figure  F-2.  Natural  sloshing  frequency  and  sloshing  modal  mass  for  lateral  sloshing  of  a 

(20'h  xlO'wx  201)  rectangular  tank. 
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APPENDIX  G 


STRAIN  AND  VELOCITY  HISTORIES  (MEASURED  AND  CALCULATED)  FOR 
THE  "RED  SNAPPER"  MODEL  TESTS 
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Figure  G-l .  Starboard-side  fiber  strain  histories  for  elements  6  and  9  for  Test  2  of 

the  "Red  Snapper"  model  test  (horizontal  and  vertical  mode  1  whipping),  and  strains 

calculated  using  the  computational  algorithm. 
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Crown  axial  strain  vs.  time,  element  6 
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Figure  G-2.  Crown  and  keel  fiber  strain  histories  for  elements  6  and  9  (respectively)  for 

Test  2  of  the  "Red  Snapper"  model  test  (horizontal  and  vertical  mode  1  whipping),  and 

strains  calculated  using  the  computational  algorithm. 
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Figure  G-3 .  Horizontal  velocity  histories  for  nodes  1 0  and  20  for  Test  2  of 

the  "Red  Snapper"  model  test  (horizontal  and  vertical  mode  1  whipping),  and  velocities 

calculated  using  the  computational  algorithm. 
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Axial  strain  vs.  time,  element  6 
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Axial  strain  vs.  time,  element  17 
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Figure  G-4    Starboard-side  fiber  strain  histories  for  elements  6  and  1 7  for  Test  3  of 
the  "Red  Snapper"  model  test  (horizontal  mode  2  whipping),  and  strains  calculated  using 

the  computational  algorithm. 
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Horizontal  velocity  vs.  time,  node  10 
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Horizontal  velocity  vs.  time,  node  20 
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Figure  G-5.  Horizontal  velocity  histories  for  nodes  10  and  20  for  Test  3  of 
the  "Red  Snapper"  model  test  (horizontal  mode  2  whipping),  and  velocities  calculated 

using  the  computational  algorithm. 
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